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EXECUTIVE  SUMMARY 


In  this  investigation  the  fields  due  to  a  current  source  embedded  inside  a  layered  media  terminated  by  a 
perfectly  conducting  ground  plane  was  studied  by  evaluating  the  Sommerfeld  integrals  along  the  real  axis. 
The  method  consists  of  locating  the  proper  surface  wave  poles  and  then  calculating  their  corresponding 
residues  for  subtracting  the  singular  behavior  of  the  integrand  at  the  pole  locations.  Furthermore,  a  new  and 
improved  closed-form  expression  to  calculate  the  rapidly  oscillatory  tail  part  of  the  Sommerfeld  integral 
was  derived  that  is  specifically  more  appropriafe  for  piecewise  mulfilayer  configurations.  The  numerical 
resulfs  were  included  for  a  single  layer  subsfrale  fhaf  showed  marked  discontinuous  behavior  in  fhe 
componenf  of  fhe  dyadic  Green’s  function  near  fhe  air-fo-subsfrafe  inferface.  Resulfs  for  laferal  separations 
p  =  lOX,  lOOX  are  included  here.  An  improved  mefhod  for  defermining  fhe  inifial  locafion  of  fhe  poles  for 
an  elecfrically  fhick  layered  media  has  been  developed.  The  final  locafion  of  fhese  poles  were  refined  by  fhe 
well-known  Newfon-Raphson  mefhod.  The  Sommerfeld  infegral  evaluafion  was  carried  ouf  by  assuming 
fhaf  fhe  infegrand  has  simple  firsf  order  poles.  If  was  however  concluded  fhaf  for  slacked  piecewise  conslanl 
layers  fhe  numerical  delerminalion  of  fhe  pole  localions  become  numerically  difficull  as  fhe  number  of  layers 
increases. 

For  a  piecewise  conslanl  multilayered  media  fhe  Green’s  function  can  assume  differenl  forms  as  fhe 
observation  poinl  moves  from  one  layer  lo  fhe  olher.  This  poses  difficully  for  real-axis  inlegralion  because 
ils  analylic  behavior  needs  lo  be  examined  for  “smoolhing”  ouf  fhe  integrand  lo  facililale  numerical  calcula¬ 
tion  of  fhe  Sommerfeld  infegral  fail.  The  procedure  Ihus  becomes  very  complex  if  a  continuously  slralified 
media  is  replaced  by  piecewise  conslanl  discrefe  layered  media.  An  allernalive  approach,  fhe  Phase-Inlegral 
mefhod  which  is  based  on  fhe  Wenlzel-Kramers-Brillouin  (WKB)  mefhod,  was  investigated  and  if  fumed  ouf 
fhaf  ils  applicalion  makes  if  feasible  lo  calculate  fhe  Green’s  funclion  for  bolh  piecewise  conslanl  or  conlin- 
uously  slralified  media.  In  addition,  fhe  Phase-Inlegral  melhod  is  valid  near  “lurning”  poinls  or  equivalenlly 
Ihe  Stakes  lines  which  are  identified  as  regions  of  difficully  in  numerical  calculations.  For  accurate  resulls 
via  Ihe  Phase-Integral  melhod  bolh  dominanl  and  subdominanl  terms  need  lo  be  included.  The  subdominanl 
term  allains  much  numerical  significance  near  Stakes  lines  while  Ihe  dominanl  term  becomes  insignificanl 
in  Ihe  close  vicinity  of  Ihose  regions.  New  asymplolic  techniques  which  can  provide  far  superior  numerical 
resulls  have  also  been  suggested  for  fulure  work  in  relation  lo  applicalion  of  Phase-Integral  melhods  lo  lay¬ 
ered  media  problems.  An  alternate  melhod,  Ihe  hybrid  ray-mode  technique,  appears  lo  have  much  potential 
specifically  for  problems  related  lo  antennas  in  embedded  in  layered  media.  In  Ibis  melhod  Ihe  Green’s 
funclion  is  decomposed  parity  inlo  modes  and  rays  which  are  mulually  related  via  Fourier  Iransforms.  The 
description  in  terms  of  bolh  rays  and  modes  caplures  Ihe  far  and  near  fields  of  radiating  sources,  and  hence 
can  provide  an  accurate  description  of  Ihe  radiation  behavior.  For  continuously  slralified  media,  Ihe  hybrid 
ray-mode  Green’s  function  can  incorporate  Ihe  Phase-Integral  solution  and  Ihus  remains  valid  near  Stakes 
regions. 

The  investigations  and  an  extensive  survey  of  Ihe  lileralure  suggesls  lhal  various  olher  techniques  such 
as  real-axis  integration,  discrete  complex  image  melhod  (DCIM),  finite-difference  lime-domain  (FDTD), 
and  a  hosl  of  asymplolic  techniques  are  necessary  lo  solve  Ihe  problem.  The  DCIM  and  FDTD  have  Ihe 
virlue  of  generating  reference  solutions  for  appropriate  cases  of  interest  The  asymplolic  techniques  are 
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useful  but  care  must  be  exercised  in  their  mathematical  development.  In  this  context  it  was  found  that  the 
existing  asymptotic  solutions  for  Sommerfeld  integrals  are  limited  in  applicability  because  of  the  manner  in 
which  they  were  derived  in  the  first  place.  Improved  methods  for  obtaining  more  general  solutions  have  been 
identified.  Finally  some  recommendations  have  been  included  in  this  report  to  justify  further  investigations 
in  this  area. 
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STUDY  OF  SOMMERFELD  AND  PHASE-INTEGRAL  APPROACHES  FOR  GREEN’S 
FUNCTIONS  FOR  PEC-TERMINATED  INHOMOGENEOUS  MEDIA 

1.  INTRODUCTION 

Electromagnetic  and  acoustic  wave  propagation  in  layered  media,  an  intensely  studied  topic  [1-7],  still 
finds  application  in  many  problems  such  as  microstrip  antennas  [8-13].  A  survey  of  the  relevant  research 
work  shows  that  the  general  aspects  of  modeling  wave  propagation  in  layered  media  has  been  treated  by 
the  Sommerfeld  approach  [11-13]  for  antenna  problems,  or  via  the  Phase-Integral  [14-24]  and  Wentzel- 
Kramers-Brillouin  (WKB)  techniques  for  microwave  propagation  in  ionospheric  and  inhomogeneous  media 
[25].  Approaches  similar  to  the  Phase-Integral  and  WKB  techniques  have  been  studied  for  determining  the 
resonant  characteristics  of  seismic  wave  propagation  [26-28].  The  subject  of  this  investigation  is  to  examine 
the  various  methods  for  determining  electromagnetic  wave  propagation  in  layered  media  with  a  particular 
emphasis  on  Sommerfeld  and  Phase-Integral  approaches. 

The  problem  of  interest  is  shown  in  Fig.  1.  The  situation  depicted  here  could  mimic  the  case  of  em¬ 
bedded  antennas  on  tactical  platforms  that  are  physically  concealed  by  reducing  physical  visibility.  Since 
such  embedded  antennas  radiate  from  within  composites  or  other  materials  the  radiation  mechanism  can  be 
understood  by  a  full-wave  analysis  including  the  effects  of  the  layered  structure.  Antenna  arrays  #1  and  #2 
contain  aperture  elements  that  radiate  from  within  a  continuous  vertically  stratified  media  media  backed  by 
a  perfecf  elecfric  conducfor.  As  a  resulf  fhe  infrinsic  propagafion  consfanf  changes  confinuously  for  z  <  H 
and  has  a  consfanf  value  of  ko  for  z  >  H.  This  abrupf  change  in  fhe  propagafion  consfanf  causes  fhe  wave 
behavior  fo  change  abrupfly  af  fhe  heighf  z  =  H. 

From  an  analysis  viewpoinf,  layered  media  propagafion  has  been  freafed  adequafely  for  single  layer  mi- 
crosfrip  anfenna  problems  [11-13]  where  one  media  is  a  lossy,  non- magnetic  dielecfric  and  fhe  ofher  is  free 
space.  The  full-wave  solution,  based  on  fhe  Green’s  function  mefhod,  shows  fhaf  fhe  air-dielecfric  inlerface 
acfs  a  guiding  sfrucfure  for  surface  waves.  The  complefe  solufion  is  described  as  inverse  Fourier  Iransforms, 
also  known  as  Sommerfeld  infegrals.  However,  fhe  dielecfric  media  was  assumed  fo  have  consfanf  permif- 
fivify  and  nof  confinuously  sfrafified  as  in  fhe  problem  of  inferesf.  Calculation  of  fhe  Sommerfeld  integrals 
is  a  cenfury-old  problem  and  research  work  is  sfill  continuing  info  ifs  evaluation  wifh  considerafion  fo  fhis 
particular  physical  problem. 

The  mosf  significanf  difficully  wifh  calculafion  of  fhe  Sommerfeld  infegrals  is  fhaf  in  fhe  neighborhood 
of  fhe  plane  confaining  fhe  source,  fhe  infegrand  in  fhe  Sommerfeld  integral  is  weakly  convergenf  and 
oscillafory.  The  oscillafions  of  fhe  integrand  is  pronounced  as  fhe  laferal  separafion  befween  fhe  source 
and  observafion  poinf  increases.  In  addifion,  fhe  locafion  in  fhe  complex  specfral  plane  of  fhe  poles  of  fhis 
infegrand  complicafe  fhe  numerical  calculation.  For  a  single  layer  problem,  as  in  microsfrip  anfennas,  a  very 
good  discussion  on  such  issues  can  be  found  in  Ref.  [12]. 

Figure  2  describes  an  “equivalenf”  problem  for  Fig.  1.  In  fhis  sifuafion  fhe  radiafion  from  any  one  of 
fhe  arrays  in  Fig.  1  is  modeled  by  fhe  equivalenf  currenfs  placed  on  a  box  wifh  5  sides.  The  confinuously 
Manuscript  approved  December  14,  2009. 
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Fig.  1  —  Two  electrically  large  arrays  #1  and  #2  radiating  from  an  infinite  ground  plane  in  presence 
of  a  continuous  vertically  stratified  media  of  thickness  H,  permittivity  e(z),  and  permeability  |J,(z). 
The  characteristic  wavenumber  of  the  continuously  stratified  medium  is  k(z)  =  CO  Aye(z)|J,(z).  Note 
that  the  intrinsic  wavenumber  is  ko  =  CO^Eolio  for  z  >  H,  where  £0  =  8.854185  X  10  Farads/meter 
and  po  =  1.256  x  10^*^  Henries/meter  are  the  constitutive  parameters  of  free-space. 
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Fig.  2  —  Representation  of  the  continuously  stratified  media  by  discrete  slabs  of  finite  thicknesses. 
The  radiating  array  has  been  replaced  by  a  equivalent  region  of  five  surfaces  with  radiating  currents 
on  the  surfaces.  Portions  of  the  region  can  extend  beyond  the  layered  region. 
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Stratified  (or  layered)  media  in  Fig.  1  is  replaeed  by  multiple  slabs  of  varying  thieknesses  and  eaeh  with 
different  but  eonstant  permittivity  and  permeability.  The  observation  point,  R,  is  now  arbitrarily  loeated 
as  shown  in  Fig.  2.  Thus,  to  ealeulate  the  field  af  P  due  fo  equivalenf  eurrenfs,  one  needs  fo  find  fhe 
Green’s  funefion  for  a  mulfilayered  media  where  fhe  number  of  layers  eould  possibly  be  large  buf  finile. 
The  propagator  (ofherwise  known  as  fhe  Green’s  funefion)  for  fhis  diserefe  mulfilayer  topology  has  been 
reporfed  fo  be  a  good  approximation  to  fhe  aefual  layered  media  whieh  is  eonfinuously  sfrafified  [1]. 

Unforlunalely,  for  more  fhan  fwo  layers  fhe  mafhemalieal  form  of  fhe  Green’s  funefion  furns  ouf  ex- 
eeedingly  eomplex  [2,  11].  In  fael,  as  fhe  poinf  P  moves  from  inside  fhe  layered  media,  erossing  fhe  various 
layers  and  info  free-spaee,  fhe  Green’s  funefion  (whieh  also  eonfains  fhe  Sommerfeld  infegral)  ehanges  [2]. 
This  ehange  is  effeeled  by  fhe  ehange  in  fhe  infegrand  of  fhe  Sommerfeld  infegral,  and  henee  fhe  Green’s 
funefion.  Thus  one  of  fhe  primary  diffieulfies  in  ealeulafing  fhe  mulfilayered  Green’s  funefion  is  fhe  ne- 
eessify  of  aeeurafely  ealeulafing  fhe  mulfiple  Sommerfeld  infegrals  fhaf  resulf  when  one  approximafes  a 
eonfinuously  sfrafified  media  by  diserefe  layers. 

The  solufion  fo  fhe  Green’s  funefion  for  a  eonfinuously  sfrafified  media  ean  be  obfained  by  fhe  WKB  [1- 
3]  or  equivalenfly  Phase-Infegral  [16-18]  mefhods.  The  Phase-Integral  mefhod  solves  fhe  one-dimensional 
wave  equation,  represenfafive  of  fhe  direefion  of  sfrafifiealion,  by  eonsidering  “fuming  poinfs.”  Physieally, 
fhe  one-dimensional  wave  equation  is  equivalenf  fo  fhaf  of  a  one-dimensional,  non-uniform  transmission 
line.  Methods  of  obtaining  the  solutions  have  been  rigorously  developed  in  Refs.  [16-20].  Compared  to 
the  diserete  multilayer  approximation  of  the  eontinuously  stratified  media,  fhe  Green’s  funefion  direefly  ob¬ 
fained  for  a  eonfinuously  sfrafified  media  does  nof  assume  various  differenf  forms.  However,  fo  aeeounf 
for  eonfinuous  uniaxial  slralifealion,  speeial  non- elementary  funelions  [17]  appear  inside  infegrand  of  fhe 
Green’s  funefion  (fhe  Airy  funefions  fhemselves  exhibif  eomplefely  differenf  behavior  for  posifive  and  nega¬ 
tive  real  argumenfs  and  henee  appear  appropriafe  fo  deseribe  wave  propagafion  in  sueh  sifuafions).  Depend¬ 
ing  on  fhe  nafure  of  fhe  sfrafifiealion,  fhe  Phase-Infegral  solufion  [16]  may  eonfain  ofher  non-elemenfary 
forms  sueh  as  Weber’s  parabolie  eylinder  or  Whiffaker’s  eonfluenl  hypergeomefrie  funelions  [18,  20]. 

From  a  eompulalional  viewpoinl  handling  one  single  Green’s  funefion  aeeounling  for  eonfinuous  slral- 
ifiealion  is  nof  mueh  more  eompliealed  fhan  evaluating  several  Green’s  funelions  for  diserete  mulfilayered 
media.  This  perspeelive  suggesls  fhaf  a  eomparison  belween  fhe  fwo  approaehes  be  done  for  an  appraisal  of 
fhe  various  approaehes.  To  fhaf  end,  fhe  primary  purpose  of  fhis  investigation  is  to  eompare  Green’s  fune- 
lions  eonlaining  fhe  Sommerfeld  infegral  for  diserefe  mulfilayered  media  versus  fhe  aefual  Green’s  funefion 
for  a  eonfinuously  sfrafified  media. 

The  earliesl  appliealion  of  Sommerfeld  Iheory  to  submarine  eommuniealion  was  reporfed  in  Ref.  [29]. 
The  radialion  behavior  of  a  dipole  immersed  inside  sea-wafer  (a  eondueling  medium)  was  explained  by  fhe 
behavior  of  fhe  Sommerfeld  infegral  aboul  fhe  braneh-eul.  This  problem  was  espeeially  dislinel  from  fhe 
Iradilional  Sommerfeld  solufion  [30].  Compulalional  slralegies  for  eleelromagnelie  wave  propagafion  in 
mulfilayered  media  and  slruelures  depends  heavily  on  analylie  approaehes  [31-32].  In  fael,  fhe  mefhod  of 
momenls  formulation  as  oullined  in  Refs.  [33,  34]  emphasizes  fhe  imporlanee  of  fhe  analylie  form  of  fhe 
mulfilayered  Green’s  funefion,  beeause  analylie  approaehes  provide  subslanlial  savings  in  eompulalional 
resourees.  In  Ref.  [35]  fhe  analylie  approaeh  fo  Green’s  funefion  formulalion  for  various  eleelromagnelie 
boundary  value  problems  is  available.  The  foeus  of  Ref.  [35]  is  fo  develop  a  eleelrie  nelwork  approaeh 
in  deriving  fhe  Green’s  funelions.  This  is  aeeomplished  by  synlhesizing  fhe  Ihree-dimensional  Green’s 
funefion  from  ils  Ihree  one-dimensional  Green’s  funelions.  These  one-dimensional  Green’s  funelions  are 
fhe  solutions  to  fhe  one-dimensional  Slurm-Liouville  problems.  In  Ref.  [35]  fhe  one-dimensional  Green’s 
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functions  are  interpreted  as  solutions  to  the  transmission  line  equations  [36,  chap.  7;  37,  chap.  8].  The 
final  methodology  involves  use  of  a  double  contour  integral  [35,  p.  285,  Eq.  37;  38,  p.  86,  Eq.  76].  This 
“network”  approach  is  particularly  helpful  for  stratified  or  piecewise  consfanf  media,  and  is  the  preferred 
approach  studied  here  in  Chapter  3. 

Some  relevant  work  involving  the  Green’s  function  for  layered  media  [31,  32]  involves  asymptotic 
evaluation  of  integrals  of  the  various  inverse-Eourier  transforms.  Such  methods  traditionally  involve  rigorous 
application  of  analytic  function  theory  to  boundary  value  problems.  General  methods  of  treating  various 
boundary  value  problems  in  acoustic  [4]  and  oceanic  [5]  media  can  generally  be  found  in  Refs.  [39-46] 
which  also  includes  various  useful  practical  techniques  of  asymptotic  evaluation  of  integrals  [42,  43]  for 
numerical  calculation  of  the  Green’s  functions. 

The  study  of  efficient  numerical  approaches  for  calculation  of  the  Green’s  function  for  microstrip  and 
other  layered  media  has  been  studied  extensively.  The  selected  Refs.  [47-59]  contain  information  of  con¬ 
tinuing  importance.  The  information  gleaned  from  these  sources  indicate  considerable  effort  in  the  direct 
numerical  calculation  of  the  Sommerfeld  integral  without  resorting  to  approximate  (asymptotic)  techniques 
[31,  32].  In  particular,  the  real  axis  integration  of  Sommerfeld  integrals  [56,  57]  is  extensively  employed  in 
Chapter  2  for  the  numerical  calculation  of  the  Green’s  function  for  the  vertical  electric  field  of  a  horizonfal 
elecfric  dipole  (G^x)-  The  difficully  of  real  axis  infegration  is  realized  when  fhe  relafive  verfical  separafion 
befween  source  and  observer  locafions  decreases.  This  difficully  is  furlher  compounded  by  the  fact  that  for 
mutual  coupling  between  elements  in  a  microstrip  antenna  array,  the  lateral  separation  could  be  very  large 
and  oscillatory  behavior  of  the  integrand  in  the  Sommerfeld  integral  causes  divergence.  To  that  end  Ref.  [51] 
has  developed  some  closed  form  approaches  for  calculation  of  the  “infinite”  tail  part  of  the  Sommerfeld  inte¬ 
gral  but  the  approach  applies  to  piecewise  constant  media  and  is  limited  by  the  fact  that  when  the  observation 
point  moves  vertically  from  one  layer  to  the  other  (see  Eig.  2)  the  changes  in  electrical  characteristics  due 
the  different  individual  discrete  layers  is  not  properly  captured.  Various  closed  form  results  involving  Bessel 
functions  [60]  with  judicious  choice  of  numerical  integration  [61]  and  formulas  for  Eaplace  transforms  [62] 
appear  useful  in  numerical  calculation  for  multilayer  problems. 

In  contrast  to  the  direct  numerical  calculations,  such  as  in  Ref.  [57],  computationally  efficient  forms  for 
Sommerfeld  integrals  for  specialized  source  and  observer  locations  have  been  developed  for  PEC -backed 
single  and  double  layer  dielectric  substrates  [63-65].  Though  these  are  of  limited  utility,  they  are  compu¬ 
tationally  extremely  efficient  within  their  domains  of  validity.  In  particular  the  result  in  Ref.  [63]  has  been 
reported  to  reduce  the  computational  time  of  a  full-wave  Method  of  Moments  (MoM)  formulation  by  several 
orders  of  magnitude  when  applied  to  electrically  large  microstrip  arrays  that  have  very  widely  separated  ele¬ 
ments.  However,  these  solutions  cannot  be  generalized  to  arbitrary  topologies  as  was  observed  in  Refs.  [31, 
32].  The  other  limitation  of  these  asymptotic  solutions  is  that  they  do  not  apply  to  multilayer  media  having 
more  than  two  layers  with  arbitrary  electrical  constitutive  parameters.  In  fact,  the  double  layer  solution  in 
Ref.  [64]  has  restrictions  on  the  relative  vertical  separations  between  source  and  observer  locations.  Regard¬ 
less  of  these  restrictions,  as  noted  in  Ref.  [34],  unless  some  form  of  analytic  reduction  is  worked  out  for  the 
Sommerfeld  integrals  a  MoM  strategy  becomes  computationally  expensive  particularly  at  high-frequencies. 

An  interesting  approach,  the  Discrete  Complex  Image  Method  (DCIM)  [66-68],  formally  avoids  the  di¬ 
rect  calculation  of  the  Sommerfeld  integrals  because  it  uses  the  Sommerfeld- Weyl  identity  and  after  some 
mathematical  manipulations  expresses  the  entire  result  as  a  finite  series  of  exponentials.  The  coefficients  of 
the  series  are  calculated  by  the  matrix  pencil  method  in  Refs.  [67,  68]  based  on  Ref.  [59].  This  is  a  particu¬ 
larly  appealing  alternative  and  appears  to  have  much  promise  in  solving  multilayer  problems.  However,  the 
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extension  of  DCIM  to  eontinuously  stratified  media  does  not  appear  to  be  available.  The  DCIM  approaeh 
ean  furnish  the  referenee  solution  whieh  ean  be  eompared  to  solutions  from  other  analytieally-dominant 
teehniques.  A  review  of  other  reeent  approaehes  reported  in  Refs.  [69-73]  suggests  a  renewed  emphasis 
on  obtaining  elosed-form  solutions.  However,  it  is  interesting  to  note  that  elosed-form  solutions  ean  also 
be  eomplieated  in  their  mathematieal  form  and  may  require  evaluation  of  Lipsehitz-Hankel  integrals.  Also, 
Refs.  [66-73]  do  not  address  the  eontinuous  stratifieation  ease. 

For  either  real-axis  integration  or  asymptotie  evaluation  of  Sommerfeld  integrals  there  exists  poles  in 
the  integrand  whose  loeation  and  residues  need  to  be  evaluated.  The  teehnique  pursued  in  Chapter  2  utilizes 
the  ealeulation  of  the  residues  upon  loeating  the  poles  for  direet  numerieal  ealeulation  of  the  Sommerfeld 
integral  (this  teehnique  is  based  on  the  approaeh  suggested  in  Ref.  [12]).  To  that  end,  there  are  various 
eonsiderations  regarding  the  type  of  the  pole.  Most  of  the  researeh  has  eentered  around  mierostrip  antennas 
whieh  are  eleetrieally  thin  planar  topologies.  As  discussed  in  Refs.  [6,  12,  13]  there  exists  one  proper  TM 
surface  wave  pole,  because  the  thickness  of  the  PEC -backed  substrate  for  mierostrip  problems  is  around 
0.01  X.  However,  as  the  thickness  increases  the  number  of  poles  that  need  to  be  considered  increases. 
Additionally,  for  asymptotic  evaluation  of  Sommerfeld  integrals  [63-65]  the  real-axis  contour  is  deformed 
and  made  to  pass  from  the  top  (proper)  to  the  bottom  (improper)  Riemann  sheet,  because  of  the  presence 
of  a  branch  cut  at  the  free-space  wavenumber  ko.  In  doing  so  improper  poles  lying  on  the  bottom  sheet 
may  be  captured  or  crossed  by  the  steepest  descent  path.  The  residues  at  these  improper  poles  are  needed 
for  accurate  evaluation.  Thus  pole  location  plays  a  very  important  role  and  representative  techniques  are 
reported  in  Refs.  [74-79].  However,  for  electrically  thick  substrates  the  information  in  Ref.  [52]  appears 
very  relevant.  Sophisticated  techniques  have  been  reported  in  Ref.  [75]  but  as  expected  its  implementation 
can  consume  significant  computation  times. 

For  layered  media  a  computationally  efficient  form  the  Green’s  function  is  available  in  Refs.  [80-85]. 
This  technique  synthesizes  the  Green’s  function  in  terms  of  finite  number  of  rays  and  modes  through  the 
use  of  the  Poisson  sum  formula.  The  process  involves  asymptotic  evaluation  of  integrals  [42,  43].  This 
ray-mode  approach  in  deriving  the  Green’s  function  appears  physically  most  appealing  for  layered  media 
because  a  layered  media  can  have  various  types  of  modes  (trapped,  surface,  leaky,  lateral  waves)  plus  rays. 
Often  one  form  (ray  or  mode)  is  inadequate  and  a  combination  of  both  is  therefore  a  judicious  choice.  This 
ray-mode  approach  has  been  pursued  for  only  few  practical  problems  in  Ref.  [85]  and  hence  it  would  be 
worth  investigating  its  application  to  the  topology  shown  in  Fig.  1. 

Finally,  Refs.  [86-100]  provide  new  aspects  on  asymptotic  evaluation  of  integrals,  which  are  at  the  heart 
of  development  of  computationally  efficient  techniques  for  layered  media  problems.  The  information  in 
these  papers  is  more  advanced  than  that  available  in  Refs.  [42,  43].  The  reason  is  that  a  new  technique 
hyperasymptotics  has  been  developed  in  Refs.  [88,  90,  95]  whose  salient  feature  is  to  provide  numerical 
accuracy.  In  the  past  the  techniques  in  Refs.  [42,  43]  focused  on  obtaining  closed-form  solution  to  the 
various  integrals  without  considering  the  accuracy  of  the  final  resulf.  Thus,  discrepancies  could  nol  be  easily 
explained  when  numerical  comparisons  were  performed  befween  exacf  and  asympfolic  solutions.  Often  fhe 
errors  could  nol  be  accounled  for  because  of  malhemalical  inlraclabilily.  The  informalion  in  Refs.  [99,  100] 
specifically  addresses  fhis  problem.  If  is  shown  fhaf  fhe  numerical  errors  in  asympfolic  expansion  happen 
around  Slokes  lines  when  fhe  dominanl  and  subdominanl  solulions  inlerchange  roles.  These  references  are 
included  here  because  of  Iheir  polenlial  ufilily  in  fulure  invesligalions. 

The  conlenls  and  conlribufions  of  fhe  various  chaplers  in  fhis  reporl  are  briefly  discussed  nexl.  In  Chap¬ 
ter  2,  a  new  improved  melhod  for  direcl  real-axis  inlegralion  of  Sommerfeld  inlegral  fail  is  developed.  Unlike 
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Ref.  [51],  this  method  ean  be  applied  to  diserete  multilayer  media  and  this  approaeh  ean  aeeurately  eapture 
the  physieal  eharaeteristies  of  the  problem.  This  ehapter  eontains  implementation  of  the  algorithm  utilizing 
the  information  in  Ref.  [52]  to  loeate  the  proper  surfaee  wave  poles  (for  real-axis  integration,  the  path  of  in¬ 
tegration  stays  always  on  the  top  sheet  whieh  is  why  only  proper  surfaee  wave  poles  that  satisfy  Sommerfeld 
radiation  eondition  need  to  be  eonsidered). 

In  Chapter  3,  the  overall  approaeh  to  the  eontinuously  stratified  media  is  developed  following  eontour 
integral  representations  [35].  The  WKB  and  Phase-Integral  solutions  for  the  eharaeteristie  Green’s  funetion 
along  the  z-stratifieation  is  discussed  primarily  from  the  information  available  in  Ref.  [17].  It  is  shown 
through  some  simple  examples  that  the  Phase-Integral  method  is  a  general  case  when  conventional  WKB 
solutions  fail  at  the  turning  points.  The  situation  is  salvaged  by  generalizing  the  solution  to  include  non¬ 
elementary  special  functions  [16-20]  that  remain  valid  at  the  turning  points  and  reduce  to  the  conventional 
WKB  solutions  when  the  observation  point  moves  away  from  the  turning  points  (or  Stokes  lines).  The 
complete  specific  solution  is  nol  obfained  buf  fhe  overall  approach  is  discussed  here. 

Chapfer  4  discusses  fhe  various  approaches  and  confains  suggesfions  regarding  fulure  work.  This  is 
followed  by  conclusions  of  fhis  investigation  and  an  extensive  lisf  of  references  followed  by  some  imporfanf 
appendices  relevanf  fo  fhe  various  aspecfs  of  fhe  informafion  in  fhis  reporf. 


2,  Gjx  GREEN’S  EUNCTION  EOR  SINGLE  AND  DOUBLE  LAYER  PEC-TERMINATED 
MEDIA 

In  fhis  chapter  fhe  Green’s  function  for  single  and  double  layer  PEC-ferminafed  subsfrales,  shown  in 
Fig.  3,  is  calculafed  for  large  lateral  separafions.  The  calculafions  follow  fhe  fradifional  real  axis  infegrafion 
of  Sommerfeld  integrals  [11,  56-59]  as  shown  in  Fig.  4  wifhouf  any  recourse  fo  special  confour  deformations 
or  asympfolic  formulalions  [55,  63-65].  The  advanfage  of  real  axis  infegrafion  is  fhaf  fhe  process  is  fracfable. 
However,  for  exfremely  large  laferal  separafions  on  fhe  order  of  fhousands  of  wavelengfhs  fhe  existing 
mefhods  involving  real  axis  infegrafion  become  numerically  unsfable.  This  is  due  fo  fhe  divergenf  oscillafory 
behavior  of  fhe  Sommerfeld  integrand.  The  convergence  behavior  becomes  very  slow  when  fhe  observation 
poinf  is  close  fo  fhe  interfaces  in  Fig.  3. 

The  scope  of  fhe  numerical  resulfs  is  limifed  fo  fhe  calculafion  of  only  fhe  G^x  componenf  of  fhe  dyadic 
Green’s  funclion  for  single  layer  media.  Specific  numerical  resulfs  for  fhe  double  layer  media  could  nol 
be  included  due  fo  time  consfrainls.  The  resulfs  were  generated  using  fhe  soflware  MATHEMATICA  for 
fhe  purpose  of  demonslraling  fhe  feasibilily  of  real-axis  infegrafion  for  large  laferal  separafions.  In  general, 
fhe  dyadic  Green’s  funclion  due  fo  a  horizonlal  eleclric  dipole  has  nine  componenls.  For  a  mullilayer 
media,  each  of  Ihese  componenls  has  complicaled  malhemalical  forms  due  fo  fhe  presence  of  fhe  infinile 
inlegrals  involving  Bessel  function  plus  multilayer  effecls  in  fhe  inlegrand.  To  explicilly  evaluate  all  fhe  nine 
componenls  of  fhe  Green’s  funclion  is  Iherefore  a  formidable  lask  and  for  lhal  reason  allenlion  is  directed  in 
fhis  reporf  fo  illuslrale  fhe  fealures  of  fhe  general  melhod  of  compulalion  for  fhe  G^  componenf  only. 

The  primary  objective  is  fo  generale  a  reference  solulion  for  G^x  which  can  be  compared  againsl  various 
ofher  lechniques  fo  be  invesligaled  laler.  The  presenl  melhod  of  real-axis  infegrafion  for  calculating  fhe 
single  and  double  layer  G^  Green’s  funclion  componenf  involves  fhe  following  steps: 


Study  of  Sommerfeld  and  Phase-Integral  Approaches  for  Green’s  Functions 


1 


-  -  -  +  oo 

Fig.  3  —  Ax  oriented  horizontal  electric  current  element  (Hertzian  dipole)  radiating  on  the  air- 
dielectric  interface  of  a  single  layer  and  from  inside  a  double  layer  media,  both  terminated  by  a  PEC 
ground  plane. 


(1)  Divide  the  real-axis  integration  into  three  distinet  regions:  (a)  0  <  ^  <  kp,  (b)  kp  <  ^  <  and  (e) 

^  <  +°°-  Determination  of  the  “breakpoint”  ^a  [57]  is  most  eritieal  for  real  axis  integration 
(diseussed  later  in  this  ehapter). 

(2)  Loeate  the  proper  surfaee  wave  poles  that  are  elose  to  the  integration  path  on  the  real  axis,  ealeulate  the 
residues  at  these  poles,  and  subtraet  the  effeets  of  the  pole  singularities  from  the  original  Sommerfeld 
integrand.  These  proper  surfaee  wave  poles  lie  in  region  (b)  of  step  (1). 

(3)  Identify  the  asymptotie  behavior  of  the  Sommerfeld  integrand  for  the  speetral  variable  ^  ^  +oo;  it 
turns  out  that  the  asymptotie  behavior  is  appropriate  in  region  (e)  of  step  (1).  This  is  ealled  the  “tail” 
part  of  the  integration  whieh  ean  be  evaluated  analytieally.  This  further  eomprises  of  the  following 
steps,  assuming  the  determination  of  the  “breakpoint”  ^a  [57],  as  below: 

(a)  Numerieal  integration  of  the  portion  of  the  tail  from  9?e(ki )  <^< 

(b)  Analytie,  elosed-form  integration  of  the  remainder  of  the  “tail”  from  ^a  <^<  -1-00. 
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Fig.  4  —  Sommerfeld  contour  for  9le(^)  integration.  The  path  of  integration  always  stays  on  the 
top  (proper)  Riemann  sheet.  All  poles  of  the  Sommerfeld  integrand  for  G^x  lie  with  the  range: 
k()  <  ^  <  9le[koy^].  These  poles,  denoted  by  *  p  it*  the  figure,  lie  on  the  top  (proper)  sheet  and 
hence  are  proper  surface  wave  poles  satisfying  the  Sommerfeld  radiation  condition. 


These  aspeets  are  detailed  in  the  seetions  below.  The  discussion  begins  with  a  single  layer  PEC -terminated 
substrate  followed  by  the  double  layer  problem. 

2,1  Analytic  Nature  of  the  Green’s  function 

This  section  begins  with  the  formulas  for  the  Green’s  function  for  a  single  layer  followed  by  a  double 
layer  geometry.  The  objective  is  to  identify  the  analytic  nature  of  the  Sommerfeld  integrands.  No  attempt 
is  made  to  perform  contour  deformation  subsequently  followed  by  steepest-descent  methods  traditionally 
found  in  the  literature.  Instead,  our  primary  focus  is  on  real-axis  integration  of  the  Sommerfeld  integral 
[56].  The  real-axis  integration  will  have  only  proper  surface  wave  poles  lying  on  the  proper  Riemann  sheet, 
i.e.,  the  sheet  on  which  the  Sommerfeld  radiation  condition  is  satisfied.  The  leaky-wave  poles,  which  do 
not  satisfy  the  Sommerfeld  radiation  condition,  are  not  encountered  when  performing  real-axis  integration 
because  the  real-axis  integration  range,  viz-,  0  <  ^  C  o®,  always  stays  on  the  top  (proper)  Riemann  sheet. 

However,  asymptotic  formulations,  especially  like  Refs.  [63,  64],  do  consider  leaky-wave  poles  which 
lie  on  the  bottom  (or  improper  Riemann  sheet).  These  improper  leaky  wave  poles  do  not  satisfy  the  well- 
known  Sommerfeld  radiation  condition,  but  are  encountered  when  the  real-axis  is  deformed  across  branch 
cuts  so  these  poles  are  encountered.  The  steepest  descent  contour  in  Ref.  [63]  must  cross  the  branch  cut  and 
hence  has  to  traverse  from  the  top  to  the  bottom  and  back  to  the  top  sheet  in  a  two-sheeted  Riemann  surface. 
In  such  a  traversal,  the  steepest  descent  path  encroaches  improper  poles  of  the  Sommerfeld  integrand  that  lie 
on  the  bottom  sheet.  The  authors  of  Ref.  [63],  in  their  development  of  an  uniform  asymptotic  formulation, 
considered  improper  poles  that  come  close  to  the  saddle  point  for  the  special  location  when  both  the  source 
and  observer  points  are  exactly  on  the  air-dielectric  interface. 
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The  limitation  of  the  preceding  approach  and  the  ones  pursued  here  is  that  we  assume  that  the  poles  of 
the  Sommerfeld  integrand  are  simple,  first  order  type.  Majority  of  the  literature  on  pole  location  in  layered 
media  [6,  10-12,  74-79]  involving  evaluation  of  Sommerfeld  integrals  is  restricted  to  microstrip  antenna 
problems.  Recent  work  in  Ref.  [79]  is  related  to  the  development  of  an  algorithm  for  pole  locations  in 
a  planarly  layered  media  close  to  the  topology  in  Fig.  3,  but  with  no  PEC  termination.  The  information 
gleaned  from  Ref.  [79]  indicates  that  locating  poles  (of  all  orders  and  types)  in  a  layered  media  is  still  a 
subject  of  continuing  investigation. 

Normally  the  thickness  of  microstrip  antennas  is  less  than  0. 1  X  and  for  real-axis  integration  only  one 
proper  surface  wave  pole  is  adequate  in  most  cases  [55,  56].  However,  for  arbitrary  electrically  thick  layered 
media  such  simplifications  are  inadmissable.  Unfortunately,  due  to  limitations  in  available  resources,  we 
have  not  investigated  in  detail  the  nature  of  the  poles  for  single  and  double  layer  substrates  for  arbitrary 
electrical  thickness  and  constitutive  electrical  properties. 

2. 1. 1  Green ’s  function  for  a  Single  Layer  PEC-backed  Substrate 

The  Gjx  Green’s  function  for  the  single  layer  geometry  shown  in  Fig.  3  is  directly  obtainable  from 
Ref.  [10,  pp.  53-54].  For  a  horizontal  electric  dipole,  oriented  along  the  x-axis  and  located  at  the  origin  one 
can  write 

J  =  xIo  AZ  5  (x)5  (y)5  (z)  =xp^  5  (x)5  (y)5  (z) .  (1) 


The  z-components  of  the  electric  fields  generafed  by  Eq.  (1)  are  differenl  in  differenl  layer  regions  in  Fig.  3. 
In  Eq.  (1),  xp_,.  =  xIoAZ  is  fhe  elecfric  currenf  momenf.  The  fields  can  be  calculated  by  obfaining  fhe  scalar 
Green’s  function  [35,  chaps.  2,  3,  5].  The  represenfafion  of  fhe  fields  relafing  fhe  propagafor  (Green’s 
funclion)  mafrix  fo  fhe  source  Eq.  (1)  is  given  by 


F 

^xx 

■ 

G,3, 

G«  ■ 

’  Px  ’ 

F 

^yx 

=  -j^^o 

Grr 

Gk 

• 

0 

F 

^ZX 

G^ 

G., 

G,, 

0 

(2) 


where  {Exx,^yx,^vc)  is  the  elecfric  field  generafed  by  a  x-direcfed  horizonfal  elecfric  dipole.  The  z-componenf 
of  fhe  elecfric  field,  radiated  by  a  x  orienfed  horizonfal  elecfric  dipole,  for  bofh  single  and  double  layer 
PEC -backed  subsfrales  are  available  in  Ref.  [10,  pp.  53-54,  Eqs.  3.3-3.16].  From  fhese  expressions,  fhe 
corresponding  Green’s  function  can  be  easily  found  using  Eq.  (2).  For  a  single  layer  dielecfric  wifh  fhe 
horizonfal  elecfric  dipole  af  fhe  air-dielecfric  inferface  and  field  poinf  in  air  as  in  Fig.  3  fhis  reads  as  follows 


GU  =  - 


271  k; 


■cos(|)  /  ^^Ji(p^)e 


Kilan(Kid) 


S^Ko-h  jKilan(Kid) 


d^ ,  for  0  <  z  <  +oo^ 


(3) 


and  for  fhe  field  poinf  inside  fhe  subsfrafe  fhe  same  Green’s  funcfion  fakes  fhe  form 


1  -JGi 


G^  = 


271  kp 


cos(|) 


-|-oo 

f  ?'2t  'i 

cos(Ki[d- |z|])  Ko 

/  s  ■'ivPSl 

0 

COs(Kid)  SrlKo-h  jKifan(Kid) 

d^ ,  for  —  d  <  z  <  0.  (4) 
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In  Eqs.  (3)  and  (4),  the  z-directed  propagation  constants  are 

/  fork^>^, 

Km  =  ^ 

for^>k™. 

In  Eqs.  (3-5),  m  =  0  and  I  designate  the  location  of  field  points  in  air  and  inside  the  substrate,  respectively. 
We  note  that  ko = (0  ^EoP-o  is  the  instrinsic  wavenumber  in  air  andki  =koy^s)iip)T  is  the  intrinsic  wavenumber 
inside  the  substrate.  The  relative  permittivity  and  permeability  of  the  substrate  are  s^i  and  la^i ,  respectively. 
In  this  report,  and  specifically  in  fhis  chapfer,  we  have  chosen  a  non-magnefic  lossy  dielecfric  subsfrale 
(jX^i  =  1,  £ri  =  Ie^i  I  [1  —  y'fanS]).  The  loss  fangenf  of  fhe  dielecfric  is  fan5 . 

If  is  inferesfing  fo  nofe  fhaf  as  fhe  observer  poinf  P  moves  from  inside  fhe  subsfrafe  (m  =  1)  info  fhe  air 
(m  =  0),  fhe  form  of  fhe  Green’s  funcfion  changes.  The  surface  wave  poles  lie  in  fhe  range  ko  <91e(ki). 

This  aspecf  is  discussed  in  Refs.  [6,  7,  11,  38].  However,  when  ^  >  ko,  we  have  Kq  =  buf 

Ki  =  —  Subslifufing  fhese  in  fhe  denominafors  in  Eqs.  (3)  and  (4),  one  nofes  fhaf  fhe  proper  surface 

wave  poles  are  given  by  fhe  roofs  of 


Dtm(^  )  =  -£,1  y/k?-^2tan[dY^k?-^2]  =  0.  (6) 

Solufions  fo  Eq.  (6)  are  discussed  in  fhe  nexf  secfion  for  fhe  single  layer  case. 

The  infegrands  in  Eqs.  (3)  and  (4)  have  branch  poinfs  defined  by  Kq  i  =  0.  Buf  due  fo  fhe  presence  of 
fhe  facfor  Kilan(Kid)  in  bofh  numerafors  and/or  denominators  in  Eqs.  (3)  and  (4),  near  Ki  0  fhe  ferm 
Ki  fan(Kid)  dK^  =d(ki  —  fhis  indicates  fhaf  fhe  producf  is  an  even  function  of  fhe  specfral  parameter 
^  near  fhe  branch  poinf  located  af  ^  =  ±ki .  The  even  funcfion  nafure  indicafes  fhaf  fhe  near  fhe  branch  poinf 
^  =  ±ki  fhe  mulfivalued  nafure  is  non-exisfenf.  Consequenfly,  fhe  branch  poinf  af  Ki  =  0  doesn’f  confribufe 
fo  any  specific  analytic  singularify.  However,  additional  inspecfion  of  fhe  nafure  of  fhe  infegrands  in  Eqs.  (3) 
and  (4)  shows  fhaf  fhe  branch  poinf  af  ^  =  ±ko  does  confribufe,  and  is  fhe  only  confribufing  branch  poinf 
singularify.  Einally,  fhe  exisfence  of  fhe  cosine  funcfions  in  Eq.  (4)  indicates  fhaf  fhey  are  even  functions 
near  Ki  =  0  branch  poinf  locations  and  fhus  fheir  ratio  is  not  mulfivalued  near  fhe  branch  poinf  af  ^  =  ±ki . 
Nexf,  we  invesfigafe  fhe  asympfofic  behavior  of  fhe  Sommerfeld  infegrand  as  ^  ^  +oo. 

To  fhaf  end,  fhe  infegrands  in  Eqs.  (3)  and  (4)  are  idenfified  as  producfs  of  “slab”  terms  and  a  spafial 
funcfion  involving  laferal  separation  befween  observafion  and  source  location. 


cos[ki (A-\z\)]^l  (f:  .N 
cosfKid) 


poinf  P  in  Pig.  3  in  air, 

poinf  P  in  Pig.  3  inside  subsfrafe. 


(V) 


The  slab  funcfions  in  Eq.  (7)  are  defined  as 


'  Ki tan(Kid) 


Dtm(^)  ’ 

ErlKp 

Dtm(^)  ’ 


poinf  P  in  Pig.  3  in  air, 

poinf  P  in  Pig.  3  inside  subsfrafe. 


(8) 
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The  denominator,  Dtm(^  )?  in  Eq.  (8)  is  defined  in  Eq.  (6)  for  ^  >  ko,  or  ean  be  identified  easily  from  the 
term  common  to  the  denominators  in  Eqs.  (3)  and  (4). 

The  choice  of  the  breakpoint,  is  determined  by  the  asymptotic  nature  of  the  integrand  in  Eqs.  (7) 
and  (8)  that  dominates  beyond  ^  >  ^a-  In  Ref.  [57]  the  asymptotic  nature  of  the  Bessel  function  was 
utilized  for  the  weighted  average  algorithm  to  calculate  the  tail  part  beyond  the  breakpoint.  The  objective  in 
Ref.  [57]  was  to  capture  the  oscillatory  behavior  of  the  integrand  in  the  Sommerfeld  integral  (for  ^  ^  oo)  and 
thus  evaluate  the  tail  in  closed  form.  Also,  the  oscillatory  behavior  of  the  Bessel  function  was  subtracted 
out  from  the  original  integrand  and  a  subsequent  smoother  function  was  used  for  numerical  integration 
in  the  range  0  <  ^  <  ^A.  In  this  report,  the  Bessel  function  is  retained  as  it  is,  but  the  slab 

functions  are  approximated  for  ^  ^  oo.  A  closed  form  result  is  subsequently  obtained  following  some 
simple  manipulations,  using  a  variant  of  the  Sommerfeld- Weyl  identity  derived  in  Appendix  A. 

Erom  Eq.  (5)  it  follows  that  as  ^  ^  oo,  Kqj  ^  and  hence  tan(Kid)  tan(— y'^d)  =  — ytanh(^d). 
Eurther,  for  ^  ^  oo,  one  has  tanh(^d)  1.  Using  these  results  (valid  for  ^  one  has  the  following 

“asymptotic”  forms  for  the  slab  functions  in  Eq.  (8) 


>d) 


-j^x-j  ^  -j 

Erfx-y^-i-yx-y^x-y  eri+p 

ErlX-y^  _  E,l 

Erl  +  1’ 


m  =  0  point  P  in  Eig.  3  in  air, 
m  =  1  point  P  in  Pig.  3  inside  substrate. 


(9) 


Prom  Eq.  (9)  one  may  therefore  define  the  following  constants 


1 


£rl  +  1 


m  =  0  point  P  in  Pig.  3  in  air, 
m  =  1  point  P  in  Pig.  3  inside  substrate. 


(10) 


- oo 

The  interesting  aspect  is  that  in  Eq.  (10)  the  changes  discontinuously .  This  can  be  understood  if  one 
considers  for  a  pure  real  -  lossless  substrate.  As  the  observation  point  P  moves  from  inside  the  substrate 
to  air,  the  constant  in  Eq.  (10)  changes  from  pure  real  to  imaginary  for  lossless  s^i.  This  happens  rapidly 
across  the  air-substrate  interface  and  is  also  known  as  a  manifestation  of  the  Stokes  phenomenon  [88,  89]. 
Por  ^  >  ^A,  via  analytic  continuation  one  then  has 

g-iKo|z|  = 

cos[Ki(d-|z|)]  _  cosh[Y/^2T^(d-|z|)] 

~  COSh[d^^2_k2] 

In  Eq.  (11),  set  t  =  d—  |z|,  and  note  that  t  >  0.  Using  the  definitions  of  hyperbolic  functions 


cosh[tY^^2  —  kj] 


cosh[dY^^2  —  kj] 


\  |exp(t^ ^ 2  _  ^2)  ^  exp(-ty^ 

^exp(t^^2_k2)^ 
^exp(d^^2_k2); 


since  t  >  0  and  ^ 


OO 
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therefore, 


eos[Ki(d-|z|)]  ^^2) 

eos(Kid) 

Finally,  for  the  full  integrand  given  in  Eq.  (7)  one  can  “extract”  the  asymptotic  behavior  using  Eqs.  (10)  and 

(12) 

(13) 


The  locations  of  the  observation  point  P  in  air  or  inside  the  substrate  are  given  by  m  =  0  and  1,  respectively, 
as  in  Eig.  3.  The  relation  Eq.  (13)  shall  be  used  later  for  numerical  evaluation  of  the  Green’s  function. 
Next,  we  analyze  the  G^^  for  the  double  layer  case. 

2. 1.2  Green ’s  Function  for  a  Double  Layer  PEC-backed  Substrate 

The  Green’s  function  for  the  double  layer  case  with  an  horizontal  electric  dipole  can  be  written  in  the 
generic  form  as 

-j-oo 

We  also  define  the  z-directed  and  substrate  propagation  constants,  generically,  via  Eq.  (5)  with  the  index 
replacement  m^n  and  n  =  0, 1 , 2  therein  for  the  double  layer  geometry.  Here  air  {n  =  0),  layer  #  1  (n  =  1) 
and  layer  #  2  (n  =  2),  for  the  geometry  shown  in  Eig.  3.  The  integrand  in  Eq.  (14)  reads,  for  the  various 
layers,  designated  by  the  superscript,  n,  as 


[(K2d2)  +/Kodi  ) 

At  /’ 


1  i-2t  J  Sin(l^2d2) 


^vc  =  Vh{^p) 


Ai 

jcOs[K2(d2  +  z)] 
K2A1 


—  cos[Ki(di  —  z)]  +  —  sin[Ki(di  —  z) 

Zr\  Ki 


Kocos(Kidi)+  7— sin(Kidi) 


The  preceding  functional  forms  of  the  integrands  can  be  found  in  Ref.  [10,  pp. 
quantity 


for  di  <  z  C  00, 

for  0  <  z  <  di, 

for  — d2  <  z  <  0.  (15) 

53-54].  In  Eq.  (15),  the 


Ai  —  =^11  +  <^  12, 


is  given  by 


fL =  — cos(Kidi)sin(K2d2)  +  —  —  sin(Kidi)cos(K2d2), 

Zr2  £rl  ^2 

FLn  =  —  —  —  sin(Kidi)  sin(K2d2)  +  —  cos(Kidi)  cos(K2d2). 
Zr2  Kl  K2 


(16) 
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We  examine  the  braneh  point  eontributions  from  Ko,i,2  to  the  Sommerfeld  integral  Eq.  (14).  To  that  end,  one 
ean  write 

sin(K2d2)  1 

Ai  12 


where 


^11 

^12 


/  j  N  ■  /  j  N  COt(K2d2)  1 

is(Kidi)+Kisin(Kidi) - >, 

1^2  I 


(17) 


By  inspeeting  the  funetional  form  of  .^11,12  in  Eq.  (17)  one  ean  eonelude,  with  the  help  of  Eq.  (5),  that  the 
term 

sin(K2d2) 


is  an  even  funetion  of  the  speetral  parameter  ^ .  This  ean  be  inferred  by  examining  the  Taylor  series  rep¬ 
resentations  for  eos(Kidi)  and  Kisin(Kidi),  and  the  Mittag-Eeffier  form  [45,  p.  168]  for  eot(K2d2)/K2.  It 
will  turn  out  that  these  representations  are  in  even  powers  of  Ki  2  and  henee  not  multivalued  near  braneh 
points  defined  by  K1.2  =  0.  Thus,  this  term  whieh  is  also  present  as  a  faetor  in  Eq.  (15)  has  no  braneh  point 
eontributions  from  Ki  2  propagation  eonstants.  Eurther  inspeetion  of  the  funetional  behavior  of  eaeh  of  the 

integrands  as  defined  in  Eq.  (15),  show  that  the  braneh  points  at  Ky2  =  0  do  not  eontribute  to  any 

- 2 

form  of  analytie  singularity.  Eor  the  ^ ^  term  in  Eq.  (15),  it  remains  neeessary  to  examine  the  nature  {i.e., 

even  or  odd)  of  the  term  K2A1.  It  ean  be  established,  multiplying  the  various  terms  in  Eq.  (17)  by  K2  that  it 

is  indeed  an  even  funetion  of  the  term  K2.  Thus  to  all  the  integrands  in  Eq.  (15)  the  braneh  points  defined 

— 0  1  2 

by  Ki  2  =  0  do  not  confribufe.  The  only  branch  poinf  fhaf  confribufes  fo  all  ^ ^  ’  in  Eq.  (15)  is  defined  by 
K()  =  0.  Hence,  fhe  defined  by  Eq.  (14)  has  fhe  sole  branch  poinf  confribufion  from  Kq.  Since  fhe  branch 
poinf  confribufion  resulfs  in  a  confinuous  specfrum  [35,  chap.  5]  in  fhe  radiafion  zone,  fhis  implies  fhaf  fhe 
ofher  layers  acf  fo  frap  waves  af  propagafion  consfanfs  Ki  2. 

Referring  fo  fhe  lasf  equation  in  Eq.  (15),  fhe  ferm 


7COs(K2[d2  +  z]) 
K2A1 


can  be  shown  nol  fo  have  any  branch  poinf  confribufions  from  Ki  2 =0.  This  can  be  verified  by  nofing  fhaf  fhe 
denominator  K2A1  =  (,^11  +  i2)K2sin(K2d2)  is  an  even  function  of  fhe  z-direcfed  specfral  wavenumbers 
Ki,2,  because,  from  Eq.  (17)  if  follows  fhaf  ^11,12  is  an  even  function  of  Ki  2.  Thus  fhe  ratio  of  fwo  even 
funclions  is  an  even  funclion,  which  resulfs  in  fhe  deducfion  fhaf  branch  poinfs,  defined  by  Ky2  =  0  do  not 
confribufe  fo  fhe  infegrands  defined  in  Eq.  (15).  The  only  branch  poinf  confribufion  comes  from  fhe  locafion 
Ko  =  0.  The  asympfofic  behavior  of  fhe  infegrands  in  Eq.  (15)  is  described  nexf. 
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It  turns  out  that  the  quantities 

sin(K2d2)  7eos(K2[d2  +  z]) 

Ai  K2A1 

need  to  be  examined  first  for  their  asymptotie  ^  00)  behavior  sinee  these  appear  in  the  various  integrands 
defined  in  Eq.  (15).  Following  fhe  usual  feehnique  for  fhe  single  layer  ease,  one  notes  from  Eq.  (5)  fhaf 
Ki,2  ^  —7^  as  ^  ^  00.  Using  fhe  eireular  fo  hyperbolie  eosine  and  sine  fransformafions,  one  obfains  fhe 
following  sequenee  of  relafionships,  valid  for  ^  ^  00, 
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where 


C()  = 


1  +  1 

+  +  ] 

1  Gt' 

[  1 

[Erl 

Ete  J 

1  er2 

(20) 


It  is  interesting  to  note  that  the  slab  constants  are  independent  of  the  individual  physical  thicknesses  as 
^  ^  oo.  In  view  of  the  relations  Eqs.  (17-20)  one  finds  that  the  asymptotic  form  of  the  integrands  in  Eq.  (15) 
are  identical  to  that  in  Eq.  (13),  with  the  subscript/superscript  index  m  in  Eq.  (13)  replaced  by  n,  where 
n  =  0, 1 , 2  for  the  double  layer  case. 

2.2  Location  of  Poles  for  the  Integrand  of  the  Green’s  function 

In  this  section,  the  poles  of  the  Sommerfeld  integrand  for  single  layer  substrates  are  investigated.  One 
of  the  fundamental  aspects  of  real-axis  integration  of  Sommerfeld  integration  is  the  location  of  the  proper 
surface  wave  poles.  A  review  of  the  recent  investigations  reported  in  Refs.  [74-79]  suggests  the  problem 
remains  unresolved,  particularly  for  electrically  thick  substrates.  The  results  for  single  layer  are  included 
in  this  section.  Similar  analysis  for  double  layer  substrates,  though  straightforward  is  complicated  and  it 
appears  that  application  of  some  of  the  existing  methods  [79]  could  lead  to  intractability.  Thus,  the  results 
for  double  layer  are  deferred  for  future  work. 

2. 2. 1  Proper  Surface  Wave  Poles  for  Single  Layer  Substrates 

Eor  a  single  layer  substrate  the  proper  surface  wave  poles  are  given  by  an  equivalent  form  of  Eq.  (6)  that 
reads 

DTM(^)  =  -eri\/^^-koCos[(iY^ki -^2]  Y^kj -^2sin[jY/ki -^2]  =  q.  (21) 

Since  exact  analytic  solutions  to  the  above  are  not  available  [6,  11,  12,  74-79],  one  needs  to  determine  an 
‘initial  guess’  for  the  roots  and  then  use  Newton-Raphson  algorithm  (or  other  root  finding  techniques)  to 
obtain  better  estimates  for  solutions  to  Eq.  (21).  This  is  a  standard  procedure  [12].  However,  it  is  important 
to  note  that  the  convergence  to  the  desired  solution  depends  on  the  initial  guess  for  Eq.  (21). 

The  literature  contains  extensive  information  for  solutions  to  Eq.  (21)  but  for  microstrip  antennas,  which 
implies  A/X  <0.1.  Eor  layered  media,  the  preceding  condition  may  not  be  valid.  Thus  the  behavior  of  the 
integrands,  (^ ),  given  via  Eqs.  (7)  and  (8)  in  the  range  0  <  ^  C  0°  is  important.  A  practical  approach  is 
to  depict  this  information  graphically.  To  that  end,  some  representative  behavior  of  the  integrand  in  Eqs.  (3) 
and  (4)  is  presented  in  Eigs.  5  through  9  (additional  results  for  p  =  lOOOX  are  included  in  Appendix  B).  In 
all  these  figures,  depicfing  fhe  behavior  of  fhe  infegrand  in  fhe  Sommerfeld  integral,  fhe  independenf  specfral 
variable,  ^  has  fhe  unifs  of  m^^  Ofher  informafion  are  included  in  fhe  figure  capfions  fhemselves  and  hence 
are  omiffed  here  for  brevify. 

The  existence  of  fhe  proper  surface  wave  poles  is  seen  by  fhe  spikes  in  fhe  infegrand.  These  poles  occur 
in  fhe  range  ko  <  ^  <  9Ie(ki).  The  resulfs  in  Eigs.  5  fo  9  (and  Eigs.  B1  fo  B5  in  Appendix  B)  were  gener¬ 
ated  for  a  frequency  of  5  GHz.  Thus  fhe  wavelengfh  X  =  2.997925  x  10*/5.0  x  10^  =  0.0599585  meters; 
Ibis  resulfs  in  fhe  infrinsic  free-space  wavenumber  ko  =  27t/)i  =  104. 79m^^  and  9Ie(ki)  =  9Ie[ko-y/|Sri|]  = 
328.051m^^  Therefore,  all  surface  wave  poles  lie  befween  104.79  <  ^  <  328.051  for  fhe  single  layer  case 
as  shown  in  fhe  figures.  This  is  evidenced  by  fhe  spikes  in  eifher  real  or  imaginary  parfs  of  fhe  infegrand. 


16 


Chatterjee  and  Walker 


Fig.  6  —  Behavior  of  integrand  in  Eq.  (3)  for  observation  point  in  air  (region  #  0)  in  Fig.  3  for 
p  =  10?i,z  =  +?i/2,  tan8  =  0.0001,  d  =  ?i,  and  Ieh  I  =9.8 


Fig.  7  —  Behavior  of  integrand  in  Eq.  (3)  for  observation  point  in  air  (region  #  0)  in  Fig.  3.  All  data 
same  as  in  Fig.  6,  except  for  the  lossless  case  tan  8  =  0. 


Re[5"l(f)] 


Fig.  8  —  Behavior  of  integrand  in  Eq.  (4)  for  observation  point  in  substrate  (region  #  1)  in  Fig.  3  for 
p  =  10A,,z  =  -A,/2,  tanS  =0.0001,  d  =  >.,  and  IehI  =9.8 
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Fig.  9  —  Behavior  of  integrand  in  Eq.  (4)  for  observation  point  in  substrate  (region  #  1)  in  Fig.  3  for 
p  =  10>i,z  =  -?i/2,  tanS  =0.0001,  d  =  ?i,  and  Ieh I  =9.8. 
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Using  these  results  in  Eq.  (22)  we  get 


-|-00 

ErlV  Y,  -  Y 

«— 0  n—0 


v2)”+i=0. 


(24) 


Now  the  barriers  at  ^  =  ko  and  ^  =  9?e(ki)  form  the  strip  within  which  the  proper  surface  wave  poles  lie. 
This  equivalently  implies  that  the  barriers  are  defined  ‘asymptotically’  at  v  ^  0  and  u^O.  Consequently,  in 
Eq.  (24),  we  shall  have  roots/solutions  when  (r^  —  v^)  ^  0  and  v  ^  0.  This  observation  allows  application 
of  the  binomial  theorem,  viz-. 


+  00 


v(r^-vY=I  ■  (-l)"r 


p=0 
+  00 


\Pr^{n-p).,,2p+l 


(r^-v2)”+i  =  £ 


p=0 


n  + 1 


(_l)r’^2(n+l-p)^2p_ 


Using  Eq.  (25)  in  Eq.  (24)  and  truncating  both  the  infinite  series  at  the  term, 
form  of  the  equation  that  reads 


eriIv2/^+'(-l)'’£c„ 


p=0 


n—p 


N+1 


N 


-Ei-'n-iV’E 


p=0 


n—p 


n  +  1 


(25) 

we  have  the  “polynomial” 

^r2("+i-/’)=0.  (26) 


Eq.  (26)  is  solved  to  find  the  approximate  location  of  the  roots  for  v,  and  hence  the  location  of  the  p* 
pole  in  the  complex  ^  plane.  Once  is  found,  then  from  the  conditions,  viz-,  ko  <  <  9ie(ki)  and 

3m(Ko  =  Y^kg  —  ^2)  <  0,  the  appropriate  poles  can  be  chosen.  Satisfying  these  conditions  implies  the 
choice  of  poles  shall  yield  only  the  proper  surface  wave  poles,  which  is  required  for  real-axis  integration  of 
Sommerfeld  integrals. 

Once  the  determination  has  been  made  via  the  above  criterion,  the  Newton-Raphson  method  can  be 
applied  to  further  refine  fhese  pole  locations.  Eigure  10  shows  fhaf  fhe  order  of  fhe  polynomial  in  Eq.  (26) 
increases  wifh  bofh  s^i  and  elecfrical  fhickness  A/X.  The  resulf  in  Eig.  10  indicafes  fhe  difficulfy  of  defer- 
mining  relevanf  pole  locafions  for  elecfrically  fhick  subsfrafes.  Surveying  fhe  relevanf  liferafure  [74-79],  fhis 
issue  remains  an  open-problem.  The  difficulfy  is  fhaf  while  numerically  rigorous  mefhods  [75,  79]  exisf, 
fheir  implemenfafion  increases  fhe  compufafional  burden.  Erom  a  practical  perspecfive,  fhis  sfill  remains  an 
open  problem. 

The  fofal  (infeger)  number  of  proper  surface  wave  poles,  Np,  lying  in  fhe  range  ko  <  <  9?e(ki)  is 

defermined  approximately  from  fhe  relafionship 


Np  =  Int 


kody^js^il  —  1 
71 


+  1. 


(27) 


Thus,  for  Zr\  =  9.8,  d  =  X,  we  gef  Np  =  Int[5. 932958]  +  1=  5  +  1=  6  (see  Eig.  10).  The  information 
gleaned  from  fhe  resulfs  in  Eigs.  10  fo  12  may  be  summarized  fo  indicafe  fhaf  surface  wave  pole  locafions 
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Fig.  10  —  Order  N  of  the  Taylor  polynomial  equation  in  Eq.  (26)  vs.  electrical  thickness  of  the 
substrate  d/A.  for  various  substrate  permittivities. 


Fig.  11  —  Number  of  proper  surface  wave  poles  vs.  electrical  thickness  of  the  substrate  d/A,  for 

various  substrate  permittivities. 


I  ^  I  'N 
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Fig.  12  —  Number  of  proper  surface  wave  poles  vs.  substrate  permittivities  le^i  |  for  various  electrical 

thickness  of  the  substrate  d/X. 


for  arbitrarily  electrically  thick  media  is  a  difficult  computational  problem.  A  practical  guideline  needs  to 
be  determined  for  such  situations. 

The  dispersion  characteristics  of  the  proper  surface  wave  poles  with  permittivity  loss  is  shown  in  Figs.  13 
to  16.  The  total  number  of  poles  in  each  of  these  figures  remains  fhe  same,  following  Eq.  (27).  However, 
as  fhe  subsfrafe  becomes  opfically  dense,  viz-,  £,-i  =  2.22  — >  9.8,  wifh  d/X  =  0.25,  fhe  poles  can  be  seen  fo 
migrate  away  from  fhe  /kp)  axis  in  Figs.  13  fo  14. 

Similar  remarks  apply  fo  fhe  nafure  of  resulfs  Figs.  15  and  16.  Thus,  increase  of  loss  in  fhe  material 
causes  surface  wave  poles  fo  move  away  from  fhe  /kp)  axis  for  opfically  dense  or  elecfrically  fhick 
maferials.  The  movemenf  of  poles  away  from  fhe  )  axis  for  fhe  lossy  case  suggesfs  fhaf  real-axis 
numerical  infegrafion  is  much  sfable  for  fhe  lossy  case  fhan  fhe  lossless  case. 

2.3  Real  Axis  Integration  of  Sommerfeld  Integrals  for  the  Green’s  function 

In  this  section,  the  algorithm  for  calculating  the  Sommerfeld  integrals  in  the  function  in  both  single 
and  double  layer  situations  is  presented  comprehensively.  This  is  suggested  by  observing  a  lot  of  common¬ 
ality  in  their  respective  functional  behavior  over  the  range  of  integration  0  <  ^  C  oo.  There  are  three  regions 
as  identified  in  the  begining  of  this  chapter.  To  that  end,  the  two  categories  of  singular  behavior  of  the 
Sommerfeld  integrand  are:  (a)  simple  pole  singularities  and  (b)  oscillatory  behavior  over  the  tail  part  of  the 
integrand. 

We  just  rewrite  Eq.  (13),  with  a  slight  change  in  notation  as  follows 


(28) 
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Im  [i/ko] 
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Fig.  13  —  TM  surface  pole  dispersion  diagram  vs.  substrate  loss  tangents  for  a  single  layer  PEC- 
backed  substrate  with  thickness  d/X  =  0.25  and  |e|  =  2.22;  tan5  =  0.000(A),  tanS  =  0.0001(*), 
tan5  =  0.001(O)>  and  tan5  =  O.Ol(n). 


Im  Wh\ 

-0.005  : 
-0.010  : 
-0.015  : 
-0.020  : 
-0.025  : 
-0.030  - 


A 

O  1-5 


. k  Re[^//:o] 

2.0  2.5  3).0 


□ 


Fig.  14  —  TM  surface  pole  dispersion  diagram  vs.  substrate  loss  tangents  for  a  single  layer  PEC- 
backed  substrate  with  thickness  d/A.  =  0.25  and  |e|  =  9.8;  tan5  =  0. 000(A),  tan6  =  0.0001(*), 
tan6  =  0.001  (0)>  and  tan 5  =0.01(0). 
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Fig.  15  —  TM  surface  pole  dispersion  diagram  vs.  substrate  loss  tangents  for  a  single  layer  PEC- 
backed  substrate  with  thickness  d/A.  =  1.0  and  |e|  =  2.22;  tan5  =  0.000(A),  tan6  =  0.0001(*), 
tan5  =  0.001(O)>  and  tan5  =  O.Ol(n). 
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Fig.  16  —  TM  surface  pole  dispersion  diagram  vs.  substrate  loss  tangents  for  a  single  layer  PEC- 
hacked  substrate  with  thickness  d/A  =  1.0  and  |e|  =  9.8;  tan5  =  0.000(A),  tanS  =  0.0001(*), 
tan6  =  0.001  (0)>  and  tan 5  =0.01(0). 
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- oo 

In  Eq.  (28),  n  =  0, 1  for  a  single  and  n  =  0, 1,2  double  layer  substrates,  and  are  given  by  Eq.  (10)  for 
single  and  Eq.  (20)  for  double  layer  eases,  respeetively.  Next,  one  needs  to  define  the  “breakpoint”  for 
applying  the  result  in  Eq.  (28)  for  numerieal  ealeulations.  While  several  ehoiees  might  exist  [57],  for  the 
purposes  of  this  report  it  is  eonvenient  to  first  define  fhe  range:  kmin  <  ^  <  kmax  This  is  fhe  region  wifhin 
whieh  all  poles  defined  by  Dtm(^  )  =  0  (single  layer)  and  Ai  (^ )  =0  (double  layer)  eases  lie  (one  needs  fo 
refer  fo  Eqs.  (6)  and  (16)  for  explieif  funelional  forms  in  defermining  fhe  poles).  To  fhaf  end,  one  defines 


(a)  kmin  =  min  C  9ie(k„)  and  k^ax  =  max  C  91e(k„),  where,  n  =  0, 1  (single  layer)  and  =  0,1,2  (double 
layer),  as  appropriafe, 

(b)  ~  1  Akmax  =  is  chosen  as  fhe  fixed  breakpoinf  for  all  layers  instead  of  being  variable. 


Inspeefing  fhe  exaef  forms  of  fhe  infegrands  in  Eqs.  (7)  and  (15),  if  furns  ouf  fhaf  one  ean  generically  write 
bofh  single  and  double  layer  eases  as  follows 


m) 


(29) 


In  Eq.  (29),  =Dtm(^)  for  single  and  =Ai(^)  for  double  layer  eases,  respeefively.  As  slated 

before,  one  of  fhe  main  assumptions  is  fhaf  fhe  proper  surfaee  wave  poles  lying  in  fhe  region  defined  by 
kmin  <  ^  <  kmax,  are  simple,  first  order  poles.  The  residue  af  fhe  p*^  firsl  order  pole  is  given  by  [45,  p.  145] 


... 


do' 

1 

W? 

II 

iJLf 

(30) 


In  Eq.  (30) 


d^ 


jk 


eos(d^k2-^2) 


(31) 


for  fhe  single  layer  ease.  The  in  Eq.  (30)  ean  be  found  via  inspeelion  from  Eqs.  (9)  and  (15).  If 

we  assume  furlher  fhaf  Ihere  are  n  finite  number  of  sueh  poles,  Ihen  one  ean  sublrael  fhe  effeels  of  fhe  pole 
singularities  in  fhe  appropriate  region  [6,  pp.  161-164;  10;  12,  pp.  167-169].  This  requires  ealeulalion  of  fhe 
residues  al  fhe  simple  poles  via  Eqs.  (30)  and  (31).  This  ereales  a  smoolh  funelional  behavior  over  fhaf  re¬ 
gion,  and  fhe  behavior  of  fhe  Ihus  modified  inlegrand  is  shown  here.  In  addition,  beyond  fhe  breakpoinf,  one 
ean  sublrael  fhe  form  Eq.  (29)  fo  reduee  fhe  effeels  of  fhe  large  oseillalions  due  fo  fhe  Bessel  funelion.  The 
oseillalory  behavior  is  very  erilieal  when  fhe  observalion  poinl  is  near  fhe  inlerfaees,  i.e.,  z  0.  Considering 
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all  the  aspects,  we  can  write  the  integral 


Gl,  =  C'^ 


=  C”< 


'  +“ 

■0  5a  ^ 


-jKr 


^  c”  I y  +  p [3  +  37k„r  -  (V)^]  I . 


(32) 


The  last  line  in  Eq.  (32)  follows  from  Eq.  (A  10)  in  Appendix  A.  The  c”  term  is  a  constant.  The  finite  integral 
in  Eq.  (32)  expands  to 


^A 


+ 


p=i  ^  > 


d^ 


+  I 


^max  ^  4a 


kmm  k, 


(33) 


Eollowing  Ref.  [6,  pp.  161-164]  the  integral  in  Eq.  (33),  viz-. 


E 

>^max  ,  4, 

1  -In 

^max 

^min 

■JT^- 


In  arriving  at  the  preceding  result,  use  has  been  made  of  the  definition  ln(— 1)  =  —jn  [6,  p.  163]. 
Consequently,  based  on  the  results  from  the  analysis  of  the  behavior  of  the  integrands,  one  can  easily  use 
simple  numerical  integration  routines  because  the  singular  behavior  resulting  from  pole  and  oscillatory 
nature  of  the  integrands  have  been  eliminated.  In  fact,  this  is  the  fundamental  basis  of  a  wide  class  of 
techniques  for  semi-infinite  integration  of  functions  with  various  types  of  singularities.  Thus,  the  given 
in  Eqs.  (3),  (4),  and  (14)  as  improper  spectral  (inverse  Eourier  transform)  integrals  can  be  numerically 
integrated,  if  the  proper  pole  locations  can  accurately  be  calculated. 

Some  remarks  regarding  the  analytic  nature  of  the  integrands  needs  clarification  for  this  numerical  cal¬ 
culations.  Eirst,  use  must  be  made  of  Eq.  (5)  in  the  analytic  form  for  )  in  carrying  out  the  integration. 

This  ensures  that  the  analytic  continuation  is  preserved  that  is  consistent  with  Sommerfeld  radiation  condi¬ 
tion.  Second,  the  functional  form  of  (^ )  contains  an  exponent  term.  This  term  must  also  be  modified 
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in  appropriate  regions  by  the  important  analytic  continuation  result  in  Eq.  (5).  If  these  are  not  observed 
carefully,  numerical  errors  shall  accrue  and  incorrect  (and  sometimes  intractable)  answers  may  be  obtained. 

2. 3. 1  Numerical  Results  and  Discussion 

In  this  final  subsection,  the  behavior  of  the  singularity-extracted  and  original  integrand  appearing  in 
Eq.  (32),  viz.. 


)  =  ,  when  ^ 

=  ^^(^ )  -  >  ^a)  -  L  >  when  <  ^  <  k„^ax, 

p=\  s  Sp 

=  ^  ^a) >  when  ^  >  kr,^ax  =  ^a,  (34) 

is  shown  in  Eigs.  17  and  18  (additional  results  for  p  =  lOOOX  are  included  in  Eigs.  B6  and  B7).  The 
evident  importance  of  the  singularity  subtraction  is  seen  in  these  figures  and  thus  one  can  use  straightforward 
numerical  integration  methods  to  evaluate  the  resulting  integral(s)  appearing  in  the  in  Eq.  (32).  The 
results  of  real-axis  numerical  integration,  for  ^  is  shown  in  Eigs.  19  to  26.  In  these  figures,  all  data  is 
included  and  hence  are  omitted  here  to  avoid  tedium.  The  plots  are  for  observation  point  inside  (n  =  1)  to 
outside  (n  =  0)  the  substrate  and  depicted  here  as  a  function  of  z/X  for  p/X  =  10  and  100. 

To  investigate  the  computational  efficiency  of  the  singularity- subtraction  algorithm  in  Eq.  (32),  the  in¬ 
built  numerical  integration  routines  in  MATHEMATICA  were  used.  Of  course  the  algorithm  in  Eq.  (32)  was 


Fig.  17  —  Behavior  of  the  integrand  in  the  Sommerfeld  integral  for  the  G^{p ,  z,  (() )  Green’s  function 
in  Eqs.  (32)  and  (33).  Here  d  =  A,,  p  =  lOA,  (|)  =  0°,  z  =  —X/2,  tan8  =  0,  and  |eri|  =  9.8.  The 
observation  point  is  inside  the  substrate.  The  imaginary  part  of  the  integrand  identically  vanishes. 
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Fig.  18  —  Behavior  of  the  integrand  in  the  Sommerfeld  integral  for  the  Gzx{p ,  z,  (|) )  Green’s  function 
in  Eqs.  (32)  and  (33).  Here  d  =  A,,  p  =  lOA,  (|)  =  0°,  z  =  +X/2,  tan8  =  0,  and  le^il  =  9.8.  The 
observation  point  is  outside  the  substrate.  The  real  part  of  the  integrand  identically  vanishes. 


Re[G)J 


Fig.  19  —  Behavior  of  the  Gjj;(p,z,(|))  Green’s  function  vs.  z/A  when  observation  point  is  inside 
the  substrate.  Here  d  =  A,  p  =  lOA,  (|)  =  0°,  tan5  =  0.0001  and  le^i  |  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  Eor  Eq.  (32),  the 
breakpoint  =kmfl.i:  =  9^e(ki). 
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Fig.  22  —  Behavior  of  the  G^{p ,z,^)  Green’s  function  vs.  z/A,  when  observation  point  is  outside 
the  substrate.  Here  d  =  A.,  p  =  lOA,,  (|)  =  0°,  tan5  =  0.0001,  and  \Zr\ \  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  For  Eq.  (32),  the 
breakpoint  =kmax  =  9^e(ki)- 


Fig.  23  —  Behavior  of  the  Gj^(p,z,(|))  Green’s  function  vs.  z/A  when  observation  point  is  inside 
the  substrate.  Here  d  =  A,  p  =  lOOA,  (|)  =  0°,  tan6  =  0.0001,  and  |eri|  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  For  Eq.  (32),  the 
breakpoint  ^a  =kmax  =  9^e(ki). 
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Re[G«J 


Fig.  24  —  Behavior  of  the  G;._j(p,z,(|))  Green’s  function  vs.  z/A,  when  observation  point  is  outside 
the  substrate.  Here  d  =  A,  p  =  lOOA,  (|)  =  0°,  tanS  =  0.0001,  and  |eri|  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  For  Eq.  (32)  the 
breakpoint  =kmax  =  9^e(ki)- 


Eig.  25  —  Behavior  of  the  Gj^(p,z,(|))  Green’s  function  vs.  z/A  when  observation  point  is  inside 
the  substrate.  Here  d  =  A,  p  =  lOOA,  (|)  =  0°,  tanS  =  0.0001,  and  |eri|  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  Eor  Eq.  (32),  the 
breakpoint  =kmax  =  9^e(ki). 
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Fig.  26  —  Behavior  of  the  G;.;c(p,z,(|))  Green’s  function  vs.  z/A,  when  observation  point  is  outside 
the  substrate.  Here  d  =  A,  p  =  lOOA,  (|)  =  0°,  tan6  =  0.0001,  and  |eri|  =  9.8.  The  crosses  (x)  refer 
to  Eq.  (32)  and  the  triangles  (A)  refer  to  direct  MATHEMATICA  computations.  For  Eq.  (32),  the 
breakpoint  =  ^max  =  9ie (ki ) . 


also  implemented  in  MATHEMATICA.  The  differenee  in  exeeuting  the  two  eases  was  that  in  direet  numerieal 
integration  via  MATHEMATICA,  the  pole  loeations  needed  to  be  defined  as  input.  In  the  former  ease  this 
was  not  neeessary,  as  the  singularity  subtraetion  generates  a  smooth  behavior  of  the  integrand. 


The  results  in  these  figures  show  remarkably  good  agreemenf  befween  fhe  fwo  approaehes.  The  in- 
feresfing  fealure  in  fhese  resulfs  is  nofieeable  in  fhe  ehange  in  fhe  nafure  of  bofh  9ie[Gjy;(p,z,(|))]  and 
3m[Gjx(p  ,z,(|))]  near  z/X  0  or  af  fhe  inferfaee.  This  disfinefion  is  obvious  by  fhe  rapid  fransifion  in 
fhe  funefional  behavior,  and  ean  be  affribuled  fo  fhe  Sfokes  phenomenon  [88,  89]  (fhis  Sfokes  behavior  will 
be  invesfigafed  in  more  defails  for  fhe  double  layer  ease,  in  a  later  reporf).  Obviously,  if  fhere  were  additional 
layers  fhen  sueh  “diseonfinuifies”  (rapid  fransifions)  would  be  more  pronouneed  and  eapfuring  aeeurafely 
Ibis  behavior  would  be  diffieull  for  mulfiple-layers. 


The  eompufafional  differenees  in  fhe  fwo  approaehes  is  obvious  from  fhe  resulfs  in  Fig.  27,  from  whieh 
one  ean  easily  eonelude  fhaf  fhe  real-axis  infegrafion  employing  fhe  singularify  subfraefion  mefhod  would  be 
mueh  more  effieienf.  This  mefhod,  i.e.  Eq.  (32),  fakes  signifieanfly  less  lime  eompared  lo  direef  numerieal 
infegrafion  for  large  p  and  henee  appears  heller  suited  for  ealeulalions  involving  large  laleral  separations.  For 
direcf  numerieal  infegrafion,  fhe  Iradilional  real  axis  infegrafion  [56]  ean  be  modified  lo  move  fhe  portion  of 
fhe  eonlour  “slighlly  above”  fhe  )  axis  as  in  Ref.  [57]  lo  avoid  eompulalion  of  fhe  poles.  In  fhaf  ease, 
if  fhe  MATHEMATICA  soflware  was  used,  if  would  nol  require  any  a-priori  inpul  for  fhe  pole  loealions. 
However,  if  was  verified  fhaf  in  fhis  ease  direel  infegrafion  via  MATHEMATICA  would  slill  be  more  time 
eonsuming  fhen  fhe  algorilhm  in  Eq.  (32). 
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computational  time  (seconds) 


Fig.  27  —  Comparison  of  cpu  times  between  direct  MATHEMATICA  (□  —  □)  computation  and  the 
algorithm  based  on  Eq.  (32)(A  —  A)  with  variation  in  lateral  distance  The  data  refers  to  |er|  =9.8, 
tan6  =  0.0001,  A/X  =  1;  for  each  value  of  p,  the  was  calculated  for  (|)  =  0°,  and  63  datapoints 
in  the  range  —  X<z<  lOA,. 


Some  observations  regarding  the  real-axis  integration  of  Sommerfeld  integrals  via  equations  Eqs.  (32) 
and  (33)  appear  important.  Though  Eq.  (32)  is  stable  beeause  the  resulting  integrand  (as  shown  in  Eigs.  17 
and  18)  is  well-behaved  after  singularity  subtraetion,  the  oseillatory  behavior  for  ^p  ^  oo  eannot  be  eom- 
pletely  eliminated.  In  Ref.  [57]  the  asymptotie  behavior  of  Jo(^  p )  was  used  for  elosed  form  evaluation  of  the 
tail  beyond  the  chosen  breakpoint  however,  it  does  not  appear  that  the  method  is  well  suited  for  values 
of  lateral  separations  where  p  lOOOX.  This  is  because  the  oscillatory  behavior  dominates  for  very  small 
values  of  the  spectral  parameter  ^  and  hence  can  “corrupt”  any  direct  numerical  evaluation  of  Sommerfeld 
integral.  The  only  exceptions  appear  to  be  the  discrete  complex  image  method  [66-68]  in  conjunction  with 
the  matrix  pencil  method  [58,  59],  where  no  numerical  evaluation  of  Sommerfeld  integral  is  necessary. 


To  that  end,  it  appears  important  that  an  analytical  solution  (even  with  its  attendant  approximations) 
for  ^p  ^  oo  is  relevant.  Eimited  useful  analytical  formulations  are  available  [63,  64].  The  solution  in 
Ref.  [63]  is  valid  only  when  both  the  source  and  observer  are  exactly  on  the  air-substrate  interface  and  for 
single  layer  substrate.  Similarly  the  solution  in  Ref.  [64]  has  restrictions  in  terms  of  locations  of  source 
and  observer  points.  Thus  Refs.  [63,  64]  are  directly  inapplicable  for  arbitrary  variations  in  z/X  shown  in 
Eigs.  19  to  26.  But,  both  Refs.  [63,  64]  are  valid  for  large  p  separations.  To  obtain  a  single  solution  that 
is  simultaneously  valid  for  continuous  z/X  and  also  when  p  ^  oo  is  a  challenging  task  because  results  in 
this  chapter  show  that  an  accurate  solution  near  the  Stokes  surfaces  needs  to  account  for  arbitrary  electrical 
thickness  and  optical  densities  {i.e.,  arbitrary  values  of  8^  and  p^).  This  aspect  requires  application  of  more 
sophisticated  asymptotic  techniques  such  as  in  Refs.  [42-44,  87,  88,  90,  94-100]  for  improved  evaluation  of 
the  Sommerfeld  integrals  in  a  multilayered  media. 


Study  of  Sommerfeld  and  Phase-Integral  Approaches  for  Green’s  Functions 


33 


3,  PHASE-INTEGRAL  TECHNIQUE  EOR  THE  GREEN’S  EUNCTION  IN 
CONTINUOUSLY  STRATIEIED  PEC-TERMINATED  MEDIA 

The  Green’s  funetion  for  a  horizontally  oriented  Hertzian  eleetrie  eurrent  element  in  a  single  or  mul¬ 
tilayered  PEC-terminated  media  eontains  Sommerfeld  integrals,  as  shown  in  Chapter  2,  and  takes  different 
mathematieal  forms  depending  on  the  loeation  of  the  observation  position.  If  the  media  is  eontinuously 
stratified,  say  along  the  z-direetion,  the  representation  of  the  fields,  for  arbifrary  souree  and  observer  loea- 
fions,  fhrough  fhe  use  of  Sommerfeld  infegrals  beeomes  very  fedious.  The  reason  is  fhaf  fhe  eonfinuously 
slralified  media  fhen  needs  fo  be  diserefized  info  several  layers  eaeh  wifh  differenl  eomplex  wavenumbers. 
The  informalion  gleaned  from  fhe  analysis  in  Chapter  2  indieafes  fhaf  sueh  an  affempf  fo  obfain  eleefromag- 
nefie  fields  fhrough  diserefizafion  of  a  eonfinuously  layered  media  is  an  extensively  laborious  fask.  If  fhe 
number  of  layers  exeeed  fhree,  fhe  eompufafional  proeess  beeomes  very  eomplieafed. 

In  fhis  ehapfer,  a  mafhemalieally  rigorous  alfernafive  [16-20]  known  as  fhe  Phase-Integral  mefhod  is 
developed  fhaf  utilizes  fhe  eonfinuously  varying  nafure  of  fhe  wavenumber.  This  means  fhaf  for  a  z-slralified 
media  fhe  wavenumber  variation  along  fhe  z-direetion  is  eonfinuous  and  is  refained  wifhouf  any  diserefiza- 
fion  of  fhe  layered  media.  Earlier  work  [1-5,  14,  21-27]  on  wave  propagation  in  layered  media  had  exten¬ 
sively  ufilized  fhe  Phase-Integral  approaeh. 

In  whaf  follows,  fhe  eleefromagnefie  fields  due  fo  a  Herfzian  eleefrie  eurrenf  elemenf  is  derived  for  a 
medium  whieh  has  eontinuously  varying  eonsfilulive  parameters,  i.e.,  s(z)  and  |J,(z),  along  fhe  z-direetion. 
The  various  sealar  eomponenfs  of  fhe  eorresponding  dyadie  Green’s  funelions  for  a  eontinuously  varying 
media  are  also  obfained.  This  is  followed  by  fhe  eonfour  integral  represenfafion  of  fhe  fhree-dimensional 
Green’s  funetion  for  a  z-sfrafified  media.  The  one  dimensional  eharaeferisfie  Green’s  funelion  is  derived  for 
fhe  z-sfrafified  media,  and  fhis  derivation  ufilizes  fhe  Phase-Infegral  or  Wenfzel-Kramers-Brillouin  (WKB) 
mefhod  aeeounfing  for  fhe  eonfinuous  profile  of  fhe  relevanf  eonsfilulive  parameters.  The  malhemalieal 
aspeefs  of  fhe  general  formulation  for  slralified  media  are  presented  here  from  Ref.  [35],  and  from  Refs.  [16, 
17]  for  fhe  WKB  and  Phase-Infegral  solulions  fo  ordinary  differenlial  equafions. 


3.1  Green’s  Eunctions  Due  to  a  x-Directed  Horizontal  Electric  Dipole  in  a  Continuously  z-Stratified 
Media 


The  eleetrie  field  due  to  an  arbitrarily  oriented  eleetrie  and/or  magnetie  Hertzian  eurrent  element  is  given 
in  terms  of  eleetrie  and  magnetie  Hertz  potentials,  H^  from  Ref.  [35,  p.  573]  that  read 


E(r,r')  =  ^{V  X  V  X  -  7(0[t(z'){V  x  zH^}. 
s(z) 


(35) 


The  V  differential  operator  in  Eq.  (35)  operates  on  the  observation  (unprimed)  point  eoordinates.  The 
eonfinuous  variation  of  the  eonstitutive  parameters  along  the  z-direetion  indieates  the  eleetrieal  or  optieal 
inhomogeneity  of  the  ambient  medium  (a  similar  relation  for  the  magnetie  field  ean  be  found  in  Ref.  [35, 
p.  573,  Eq.  lb],  but  is  irrelevant  here  and  thus  omitted  for  brevity).  If  both  eleetrie  and  magnetie  eurrent 
elements  (Hertzian  dipoles)  are  eonsidered  as  sourees  for  the  eleetromagnetie  fields,  then  the  eorresponding 
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Hertz  potentials  are  given  by  Ref.  [35,  p.  573,  Eqs.  le.  Id] 


7cos(z')ne 

7(0|i(z')n^ 


•  V'  X  V'  X  J  -  {Mo  •  V'  X  J, 
{Jo  .  V'  X  z^^}  +  .  •  V'  X  V'  X  z^„}. 


(36) 


Unlike  V  in  Eq.  (35),  V'  operates  on  the  eoordinates  of  the  souree  (primed)  loeation.  The  •  indieates  a 
veetor  dot  (sealar)  produet.  Our  objeetive  is  to  show  how  the  preeeding  formulation  would  reduee  to  the 
general  expressions  for  the  eleetrie  fields  for  a  single  layered,  PEC-baeked  mierostrip  pateh  antenna  as  in 
Ref.  [65,  Eqs.  1-13].  In  Eq.  (36)  the  (^e,m  are  the  scalar  eleetrie  (e)  and  magnetie  (m)  potentials.  To  that  end, 
for  an  eleetrie  souree,  i.e.,  setting  Mo  =  0  in  Eq.  (36),  one  performs  the  following  algebraie  operations  for 
an  arbitrarily  oriented  Hertzian  eleetrie  eurrent  element,  Jo  =  J;  +  zJz,  that  are 


V'xV'xz^,  =  V'4^|-z(V',2C 


Jo.(V'xV'xzy  =J,. 

Jo  .  V'  X  z^„,  =  [J,  .  ( V'  X  z)  +  z  .  (V'  X  z)Jz]  U 
=  [J,.(V'xz)]^^ 

=  [J(*(V'f  xz)+jj*(z^  XZ.)]C,m 
=  -(J,xz).V',U. 


(37) 


In  the  preeeding  equations,  V  =  V,  +  z|^  and  V'  =  V',  +  z^.  In  Eq.  (37),  one  notes  the  applieation  of  the 
following  veetor  identities,  viz-. 


z.(V'xz)  =  V'.(zxz)  =0, 

I  •  (V',  X  z)  =  V', .  (z  X  J,)  =  -(J,  X  z)  .  V'„ 


has  been  used  to  simplify  some  of  the  intermediate  steps.  Eor  the  ease  of  a  x-direeted  Hertzian  eleetrie 
eurrent  element,  =  0  and  Jj  =  xP^.  Using  Eq.  (37)  one  finds  fhe  redueed  forms  for  fhe  Herfz  pofenfials 
for  a  x-direefed  Herfzian  eurrenf  elemenf  from  Eq.  (36)  are  given  by 


He  =  Pe 


1 

[jioe{z')]^  dx'dz'  ’ 


1  3^ 

7(0|a,(z')  dy'  ■ 


n  —  p 


(38) 
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The  scalar  version  of  Eq.  (35)  reads  in  terms  of  the  Hertz  potentials,  upon  its  reduction,  as 


^  z{z')d^Ue  .  , 

^  s(z)  dydz 

E  —  )  y2]-[ 


33^  ’ 


(39) 


In  Ref.  [35,  p.  573,  Eq.  2]  the  Green’s  function  for  a  continuous,  uniaxially  z-stratified  media,  corresponding 
to  scalar  Hertz  potentials  are  defined  as 

y%e  =  -jOi£{z’)%{r,r), 

y^i;m  =  -j(0[l{z')W„r{r,r').  (40) 


Substituting  Eq.  (38)  into  Eq.  (39)  for  Tig  and  interchanging  the  primed  (source)  and  unprimed  (observation) 
differential  operators,  one  can  obtain  the  following 


E  =p _ 

^  j(oe{z)j(oe{z')dxdz'^  ‘ 


(41) 


In  arriving  at  Eq.  (41)  use  has  been  made  of  the  artifice  d  /3x=  —d  /3x'.  One  can  obfain,  in  a  similar  manner, 
fhe  following  expressions  using  Eq.  (38)  fhaf  reads 


£(z03"n,  _  -Pg  ^ 

£{z)  dydz  7COS(z)  j(0s(z')  dxdy\dzdz'  j 


ja)|x(z') 


3  Tlffi 
dx 


£{z')d^ng 

s(z)  dxdz 


7C0s(z)7C0s(z')  dx^  \dzdz' J 


7(0[t(z') 


d  Tlffi 

dy 


-Pe 


ay2 


~Pe[~^j^m  + 


(42) 


In  arriving  af  Eq.  (42)  use  has  been  made  of  fhe  relationship,  d  /dt  =  —d  /dt',  where  fhe  label  t  x,y,z. 
Subsfifufing  Eq.  (42)  in  Eq.  (39),  one  obfains  fhe  elecfric  fields  due  fo  fhe  x-direcfed  horizonfal  elecfric 
dipole,  i.e.,  for  which  Jo  =  Jf  =  xPe.  Once  fhese  fields  are  obfained  fhe  scalar  componenfs  of  fhe  dyadic 
Green’s  funcfions  for  fhe  continuously  sfrafified  media  can  be  obfained  from  Eq.  (2)  in  Chapfer  2,  wifh 
Pg  =  p^  and  seffing  Py  =  0  idenfically,  fherein.  Addifionally  in  Eq.  (2),  in  view  of  Eq.  (39)  one  performs  fhe 
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following  replacements: 
Green’s  functions  reads 


■^yx 


Ejc,  Eyx  Ey,  and  E^  E^.  Then,  in  their  final  forms,  the  complete  scalar 


j 

kilo 

-j 

krio 

j 

krio 


i _ 

)7Coe(z')  j 

^(l„+ _ ^ _ TtV 

dxdy  y  7C0s(z)7C0s(z')  dzdz'  J 

_ ^ _ ^(v^rj. 

j(oe{z)j(oe{z')dxdz'^  ‘ 


jcos(z 


(43) 


These  scalar  Green’s  functions  are  valid  for  a  continuously  stratified  media,  as  well  as  for  multipart/multilayered 
media  which  has  piecewise  constant  permittivity  and  permeability,  e„  and  p.^,  respectivly,  for  the  n*  layer  of 
thickness  h„.  The  results  in  Eq.  (43)  bear  close  resemblances  with  the  x-,  y-  and  z-components  of  the  electric 
fields  excited  by  a  horizontal  electric  current  element  located  at  the  air-dielectric  interface  of  a  PEC -backed 
substrate  as  in  Ref.  [65] . 

It  is  noted  that  the  Green’s  functions  components  Gxx,yx,zx  in  Eq-  (43)  are  the  scalar  components  of  the 
general  dyadic  Green’s  functions  for  the  electromagnetic  problem  [34;  38,  chap.  2;  39].  In  general,  for  an 
arbitrary  source,  the  scalar  components  in  Eq.  (2)  are  nine  in  number,  and  all  the  nine  scalar  components 
need  to  be  found  for  this  general  case.  However,  the  scalar  potential  Green’s  functions  are  strictly  related 
to  an  arbitrarily  oriented  electric  or  magnetic  source.  Eor  the  electric  (or  magnetic)  Hertz  potential,  there  is 
the  corresponding  electric  plus  magnetic  potential  Green’s  function  [35,  p.  573,  Eqs.  Ic,  Id].  Thus,  for  an 
arbitrarily  oriented  electric  source  there  are  3  x  3  scalar  components  of  the  dyadic  Green’s  function.  But, 
regardless  of  the  source,  there  always  exists  only  two  types  of  potential  Green’s  functions  which  are,  of 
course,  scalar.  The  nine  scalar  components  of  the  dyadic  Green’s  functions  are  expressed  in  terms  of  the 
two  Hertz  potentials  or  equivalently  the  two  potential  Green’s  functions.  This  fact  is  to  be  borne  in  mind  to 
avoid  confusion  while  seeking  relationship  between  Eqs.  (38)  and  (2). 

What  we  next  need  to  calculate  are  the  Green’s  Eunctions  in  Eq.  (40)  and  the  scalar  potential 
functions  in  Eq.  (38).  The  outline  of  the  solution  to  ^e,m  Green’s  functions  is  described  below,  along  with 
the  corresponding  calculations  for  No  attempt  is  made  to  solve  the  electromagnetic  boundary  value 
problem  for  a  continuously  stratified  media  in  the  most  complete  fashion.  However,  the  general  outline  of 
the  overall  solution  is  described  here.  The  solution  procedure  for  the  continuously  stratified  media  problem 
is  for  fhe  major  pari  oblained  from  Ref.  [35,  chaps.  3,  5]  and  for  Ihe  minor  pari  from  Ref.  [38,  chap.  2]. 

One  of  fhe  main  reasons  for  fhe  above  is  lhal  fhe  profiles  for  s  (z)  and  p,  (z)  are  unknown.  Addilionally, 
fhe  general  solulion  approach  would  sfill  need  lo  be  validated  and  numerically  implemenled  wilh  fhe  slafe- 
of-arl  analysis  lechniques.  This  in  ilself  is  a  formidable  efforl  because  of  fhe  various  conslrainls  fhal  need  lo 
be  mel  for  fhe  praclical  ulilily  of  fhe  of  fhe  proposed  solulion  approach.  Thus,  such  specific  invesligafions 
are  besl  lefl  for  fulure  work  in  Ibis  area. 


3.2  Contour  Integral  Representations  of  3-Diniensional  Green’s  Functions  in  Eq,  (40) 

Erom  Ref.  [35,  chap.  3,  pp.  284-288,  Eq.  37]  the  electric  or  magnetic  potential  Green’s  functions 
can  be  synthesized  in  terms  of  three  one-dimensional  (transmission  line)  characteristic  Green’s  functions. 
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Gr(x,x';X,),G^’'”(y,y';X,) 


and  G^’'”(z,z';  Xj.),  as  follows 


^eAry)=^£  G‘:;'^{yy-Xy)GTy^-X)dXy.  (44) 


In  Eq.  (44),  =  k^(z)  —Xx  —  Xy  for  a  continuously  layered  media,  as  evidenced  by  the  continuous  variation 

of  the  wavenumber  k(z) .  The  contours  Cx,y  in  the  complex  Xx^y  planes  are  explicitly  defined  in  accordance 
with  the  singularities  of  the  individual  GJ“  functions  [35,  p.  284  ;  38,  chap.  2,  p.  86].  The  various  Xx,y  are 
the  spectral  wavenumbers  whose  explicit  relationship  to  the  one-dimensional  Green’s  functions  are  defined 
via  fhe  one-dimensional  Sfurm-Liouville  operator  equafion  [35,  p.  274,  Eq.  1]  as 


d 

ds 


q{s)+Xsw{s) 


gs{s,S-,Xs) 


(45) 


In  Eq.  (45)  fhe  generic  variable  s  (x,y)  as  appropriafe,  and  fhe  superscripfs  e,m  have  been  dropped.  If  is 
emphasized  fhaf  x  and  y  coordinafes  are  fransverse  fo  fhe  direction  of  sfrafificalion.  One  nofes  fhaf  fhe  con¬ 
tinuously  sfrafified  media  is  in  fhe  semi-infinile  domain  along  fhe  z-axis.  Boundary  conditions  are  imposed 
along  fhe  direction  of  sfrafificalion.  In  fhe  orfhogonal,  i.e.,  x-y  planes  fhere  are  no  boundary  conditions  fo 
be  mef,  and  fhe  solulions  fo  Eq.  (45)  should  physically  represenf  oufgoing  plane  waves.  The  determination 
of  GXyyz)  is  presented  in  a  separate  section. 

Eor  the  x-  and  y-component  of  the  characteristic  Green’s  function  appearing  in  Eq.  (44),  a  reduced  form 
of  Eq.  (45)  can  be  obtained  by  setting  p(x)  =  p{y)  =  p{s)  =  1,  w{s)  =  w(x)  =  w(y)  =  1  and  q{s)  =  q{x)  = 
q{y)  =  0.  The  explicit  forms  for  the  reduced  version  of  Eq.  (45)  are 

for  both  X-  and  y-  characteristic  green’s  functions.  Comparing  Eq.  (46)  with  Eq.  (45)  one  notes  that  the 
“unspecified”  eigenvalue  Xx^y  =  k^  or  k^  as  appropriafe.  Because  of  fhe  absence  of  any  boundary  condilions, 
fhe  specfrum  is  continuous,  viz-,  <  kx^y  <  one  can  assume  fhaf  fhis  problem  is  of  fhe  fype  SEP3 
[39,  chap.  2,  pp.  77-94],  wifh  a  limit  point  or  limit  circle  case.  This  involves  solving  fhe  corresponding 
homogeneous  equation  and  finding  fhe  basis  solulions  for  fhe  Iwo  regions  on  eilher  side  of  x  =  x'  (or  y  =  y'). 
Eor  fhe  corresponding  homogeneous  equation  in  Eq.  (46),  fhe  relevanl  imporlanl  delails  fo  oblain  fhe  gx 
Green’s  funclion  following  Ref.  [38,  chap.  2,  p.  64]  are  presenled  nexl.  The  modal  eigenfunctions  are 

a>i(x>x') 

<I>2(x  <  x')  = 


where  fhe  Wronskian  is 

ir(x')  =  2yk,AiA2, 
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and  the  one-dimensional  eharaeteristie  Green’s  funetions  in  Eq.  (44)  are 

Sj(x,x';X»  =  kj)  = 

=  (47) 

We  have  not  yet  determined  G^’'”(z,z';  Xj.)  in  Eq.  (44),  and  this  is  the  subjeet  of  the  next  seetion.  Before 
proceeding  further,  we  note  that  =  ^^(z)  —  k^  —  ky=K^  =  k^{z)  —  Moving  further  we  note  that  the 
one  dimensional  x-  and  y-characteristic  Green’s  functions  in  Eq.  (47)  are  the  indistinguishable  for  electric 
(superscript  e)  and  magnetic  (superscript  m)  cases,  because  the  corresponding  modal  eigenfunctions  are 
indistinguishable  [35,  chap.  3,  p.  251,  Eq.  40b].  The  difference  between  the  two  arises  in  the  calculation  of 
and  will  be  discussed  later. 

Continuing  with  the  completion  of  the  calculation  for  /)  in  Eq.  (44),  substitute  from  Eq.  (47)  and 
noting  that  dXx,y  =  d{X^y)  =  '^Xx^ydkx^y,  the  final  result  reads 


-|-oo 


Eollowing  Ref.  [35],  the  polar  transformations 

kx  =  ^  cos  a,  X  =  p  cos  (|) , 

^3,  =  Psina,  y  =  psin(|), 

dkxdky  =  p  dp  d<^, 

are  now  used  in  Eq.  (48)  to  obtain  the  following  form 

+00  +27C 

=  j  P<5?PG^’”*(z,z^P)  J  r/(Xexp{  — y'P p cos(a  —  (|) )}exp{+7P p ' cos(a  —  (|) ')}  (50) 

0  0 

where  0  <  p  <  +oo  and  0  <  (|)  <  271.  To  simplify  Eq.  (50),  without  any  loss  of  generality  of  the  physics 
of  the  problem,  one  may  conveniently  choose  the  horizontal  electric  dipole  at  the  origin  of  the  coordinate 
system,  so  that  p'  =  0  in  Eq.  (50).  In  Eq.  (50)  it  is  emphasized  that 

G,^’-(z,z';P)  ^G,^’“(z,z';XJ  =  yjk\z)-^^). 

This  means  that  G^’™  is  an  even  function  of  p.  Thus  G^’'”(z,z';— P)  =  G^’'”(z,z'; +P ).  Substituting  the 
Bessel  function  identity 

-(-oo 

exp[-yPpcos(a-^)]=  £ 


x'  =  p'cos(|)', 
y'  =  p'cos(|)', 

(49) 


(51) 
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in  Eq.  (50)  with  p'  =  0  therein,  one  ean  obtain  the  simplified  result  as  shown  in  the  following 


-|-oo  +271 

‘^eAr,r)  =  ^2  I  (7r”^^^'”^/Gr(z,z';P)J„(Pp)NP  / 


(52) 


With  the  following  result,  viz-. 


In 


271 ,  when  n  =  0 
0,  when  n  /  0 


(53) 


substituted  in  Eq.  (52),  only  the  n  =  0  term  survives  in  the  series,  and  one  finally  obfains 


=  ^  f  Gr(z,z';P)Jo(Pp)NP- 


(54) 


Eurfhermore  we  make  use  of  fhe  following  well-known  relafions  for  fhe  eylindrieal  Bessel  and  Hankel 
funelions 


Jo(Pp)  =  ^[H<'>{|Sp)  +  H<">(|Sp)] 

H<"(-z)  =  -Hf’(z)  (55) 

in  Eq.  (54)  fo  obfain  fhe  following  form  fhaf  reads  as  follows 


%Anr)  =  ^{J  Gr(z,z';P)HP)(Pp)NP  +  /  Gr(z,z';P)HW(Pp)Pr/p  1.  (56) 


The  lasf  infegral  in  Eq.  (56)  ean  be  eonvenienfly  grouped  wifh  fhe  firsl  one  by  fhe  subsfifufion  P  ^  — P  and 
furlher  use  of  Eq.  (55).  This  final  redueed  resulf  for  Eq.  (56)  reads 


+  00 

%Ar,r)  =  ^j  Gr(z,z';P)HP)(Pp)Pr/p 

—  oo 
-|-oo 

=  ^  /  Gr(z,z';P)Hf^(Pp)Pr/p.  (57) 

oog+ 


In  Eq.  (57)  fhe  lower  limif  of  fhe  lasf  infegral  indieafes  fhaf  fhe  original  real  axis  pafh  of  infegrafion  has  been 
slighfly  deformed  fo  go  anfi-eloekwise  and  also  fhis  ehange  implies  fhaf  fhe  pafh  of  integration  is  slighfly 
below  fhe  negafive  real  axis  fo  avoid  fhe  logarifhmie  braneh  poinf  singularify  of  fhe  Hq  ^  (z)  funelion  near 
z  0. 
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The  result  Eq.  (57)  allows  one  to  ealeulate  components  of  the  dyadic  Green’s  function  in  Eq.  (43). 
However  one  still  needs  the  explicit  forms  for  potential  functions  (^e,m  and  from  Ref.  [35,  pp.  573-574, 
Eqs.  5a,  5b]  they  read. 


U  =  I 


(58) 


Thus,  this  completes  the  complete  plan  for  calculating  the  fields  due  to  a  horizontal  electric  dipole  in  a 
layered  media  given  by  Eq.  (35).  The  Phase-Integral  method  discussed  next. 

3.3  Application  of  the  Phase-Integral  Technique  to  the  Green’s  Function  in  Eq,  (57) 

The  complete  solution  to  a  continuously  layered  media,  as  provided  in  Eq.  (57),  needs  to  account  for  the 
radiation  condition  and  the  boundary  condition  at  the  PEC.  To  that  end,  these  conditions  are  included  in  the 
z-stratification  and  hence  in  the  determination  of  the  one-dimensional  characteristic  Green’s 

function.  Eollowing  Ref.  [35,  p.  574,  Eqs.  5.87a,  b],  the  differential  equation 

Gz(z,z';P)  =  -5(z-z')  (59) 


d^  ,  .  d 


needs  to  be  solved,  with  appropriate  boundary  conditions,  for  g^{z,z'',^ ).  In  Eq.  (59)  the  superscripts  have 
been  dropped,  and  also 


P(z) 


,  for  electric  type, 

,  for  magnetic  type 


(60) 


Green’s  functions.  The  single  primes  in  Eq.  (60)  are  for  first  (single)  derivatives  and  q(z)  =  K^(z)  =k^(z)  — 
in  Eq.  (59).  It  is  obvious  that  the  solution  to  the  one-dimensional  Gz{z^z!\^)  Green’s  function  can  be 
related  to  the  eigenfunction  solution  to  the  corresponding  homogeneous  differential  equation  to  Eq.  (59)  that 
reads 


dz^ 


+p(z);^+q(z) 


,^(z)=0. 


(61) 


Once  appropriate  eigenfunction  solutions  to  Eq.  (61)  are  obtained,  the  G^{z,  z!\ P )  in  Eq.  (59)  can  be  obtained 
via  the  method  in  Ref.  [38,  p.  64,  Method  II].  Thus,  for  the  rest  of  this  chapter  our  focus  is  on  obtaining 
solutions  to  Eq.  (61)  rather  than  Eq.  (59). 

At  this  stage  the  complete  solution  to  Eq.  (61),  and  hence  to  Eq.  (57)  is  not  provided,  because  many 
issues  in  developing  an  engineering  solution  are  yet  to  be  specified.  Buf,  if  is  imporfanf  fo  show  how 
fhe  confinuous  sfrafificafion  of  fhe  wavenumber  k(z)  enfers  fhe  final  calculation  for  eifher  fhe  elecfric  or 
magnetic  type  Green’s  function  in  Eq.  (57).  Wifh  fhe  use  of  fhe  solution  fo  Eq.  (59),  one  can  calculafe 
fhe  fields  locafed  arbitrarily  inside  or  oufside  fhe  layered  medium.  The  advanfage  of  fhe  proposed  Phase- 
Infegral  formulafion  is  fhaf  if  is  a  single  formulafion,  and  is  valid  af  arbifrary  locations.  Unlike  fhe  piecewise 
confinuous  discrefizafion  formulafion  as  in  Ref.  [11]  for  a  fwo/fhree  layer  media,  fhe  presenf  formulafion 
does  nof  require  numerical  evaluation  of  multiple  Sommerfeld  infegrals  for  various  combination  of  source 
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and  observer  locations.  It  can  also  be  legitimately  argued  that  the  attendant  algebraic  complexities  of  using 
the  Phase-Integral  solution/representation  of  Eq.  (57)  does  provide  substantial  savings  in  computational 
resources. 


The  solution  to  Eq.  (59)  are  outlined  from  Ref.  [17,  chap.  14]  because  this  reference  provides  a  clear 
exposition  of  the  WKB  and  Phase-Integral  method  (Eanger  transformation).  The  latter,  can  legitimately  be 
called  as  a  generalization  of  the  conventional  WKB  method,  since  the  conventional  WKB  method  fails  near 
the  turning  points  [16].  The  rigorous  mathematical  justification(s)  of  the  Phase-Integral  method  is  available 
in  Refs.  [18-20],  and  are  omitted  here  for  brevity. 


Eollowing  Ref.  [17,  p.  360],  to  obtain  a  suitable  form  for  solution,  the  following  separability  condition, 
viz., 


Jff  (z)  =  ,^(z)‘^(z) 


is  substituted  in  Eq.  (61)  that,  after  some  straightforward  algebra,  results  in  the  following  equations 


=  e(z) 


dz^ 


+  [koqi(z)+q2(z)]'^(z)  =0 


where 


qi(z)  =ko 


s(z)^(z) 


q2(z) 


1 

■s"(z)' 

3 

■s'(z)' 

2 

s(z) 

“4 

s(z) 

(62) 


(63) 


The  single  and  double  primes  in  Eq.  (63)  refer  to  first  and  second  derivatives,  respectively.  The  solution  for 
the  differential  equation  in  Eq.  (62)  when  the  free-space  wavenumber  ko  ^  oo  is  the  central  subject  of  the 
Phase-Integral  or  WKB  method.  Indeed  this  is  essentially  seeking  an  high-frequency  asymptotic  solution 
[19,  20]  to  the  Green’s  function  defined  via  Eq.  (59).  We  follow  Ref.  [17,  pp.  362-364]  fo  elaborate  fhe 
salienf  fealures  of  fhe  mefhods.  The  algebraic  defails  can  be  obfained  by  following  Ref.  [17],  and  hence  fhe 
basic  fealures  of  fhe  solution  are  ouflined  nexl. 

The  solulion  is  obfained  by  selling 


‘^(z)  =exp{ko,^(ko,z)}, 

+00  r 

^(k(),z)  =  £  (64) 

in  fhe  differenlial  equation  in  Eq.  (62).  Incidenlally  fhe  approach  here,  following  Ref.  [17],  closely  parallels 
Ref.  [2,  pp.  79-88].  The  second  equation  in  Eq.  (64)  is  very  much  reminiscenl  of  whaf  is  known  as  Luneburg- 
Kline  series  in  geomelrical  optics  [2,  p.  90,  Eqs.  4.108a,  b].  Upon  performing  fhe  usual  subsfilulions  as 
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shown  in  Ref.  [17,  p.  363],  one  obtains  the  following  results  after  some  routine  lengthy  algebra 


fo  =  ±j  J  s/qfiz)dz, 

if  qi(z)  >  0; 

=  ±7  J  \/-qi(z)*, 

if  qi(z)  <  0; 

...  Cieos[ko/Vqi(z)*]-fC2sin[ko/Vqi(z)*] 

[qi(z)]^ 

if  qi(z)  >  0; 

C3  exp[ko  /  y/-qi  {z)dz\  +  C4  exp[-ko  /  y/-qi  {z)dz\ 
[-qi(z)]^ 

ifqi(z)  <0. 

(65) 

We  note  that  ^  (z)  in  Eq.  (65)  is  the  first  order  approximation,  beeause  we  have  not  used  the  higher  or¬ 
der  terms  as  maybe  apparent  from  Eq.  (64).  The  result  for  ^  (z)  in  Eq.  (65)  is  known  as  the  elassieal 
Wentzel-Kramers-Brillouin  (WKB)  approximation  [16-18].  The  higher  order  approximations  to  ^(z)  ean 
be  obtained  following  Refs.  [18-20],  whose  algebraie  forms  are  inereasingly  eomplex  in  nature.  Normally 
these  are  not  neeessary  for  a  medium  with  moderately  or  slowly  varying  eleetrieal  parameters  [2]. 

The  WKB  solution  in  Eq.  (65)  obviously  fails  when  q^  (z)  =  0,  also  known  as  the  “turning  points”  or 
zeros.  To  eireumvent  this  problem  the  Eanger  transformation  [15;  17,  pp.  375-379]  is  used.  The  feature  of 
this  approaeh  is  that  one  single  solution  ean  be  obtained  that  would  qualitatively  mateh  the  WKB  solutions  in 
Eq.  (65)  on  either  sides  of  a  turning  point.  Rather  than  elaborate  the  method  with  mathematieal  rigor,  whieh 
is  readily  available  in  Refs.  [19,  20],  an  intutitive  physieal  approaeh  is  pursued  here,  following  Ref.  [17]. 

Examining  the  solution  in  Eq.  (65)  and  the  original  differential  equation  in  Eq.  (62),  we  observe  that 
when  qi  (z)  <  0,  ^  (z)  has  a  oseillatory  (sinusoidal)  behavior;  when  qj  (z)  >  0  the  eharaeter  of  ^  (z)  ehanges 
sueh  that  it  has  exponential  growth  plus  deeay  behavior.  Beeause  in  the  lowest  order  solution  for  ^  (z),  q2(z) 
does  not  appear  in  Eq.  (65),  so  we  ean  examine  the  qualitative  behavior  by  obtaining  a  simplified  (redueed) 
form  of  the  differential  equation  in  Eq.  (64)  by  setting  q2(z)  =  0  and  qi(z)  =  ±1.  In  that  ease  we  have  two 
differential  equations  that  are  simplified  versions  of  Eq.  (64)  and  read 

d^ 

-^-fko'^(z)  =0,  for  qi(z)  =  1; 
dz 

— ^ — ko'^(z)  =  0,  for  qi(z)  =  —  E  (66) 

dz 

The  generie  solutions  fo  fhe  fwo  differenlial  equafions  in  Eq.  (66)  fhus  read 


^(z) 


ao  eos(koz)  -|-  bo  sin(koz) ,  for  qj  (z)  =  1 
ai  exp(koz) -|-bi  exp(— koz),  forqj(z)  =  — L 


(67) 


The  qualifafive  nafure  of  fhe  solutions  in  Eq.  (67)  is  fhe  same  as  in  Eq.  (65).  By  fhe  fundamenfal  fheorem 
of  algebra,  if  a  funelion  ehanges  signs  fhen  if  musf  vanish  somewhere  in  befween  fhe  fwo  eonseeufive  sign 
ehanges,  and  Ibis  vanishing  poinf  is  fhe  zero  of  fhe  funelion.  Thus  when  q^  (z)  =  -|-1  ^  —  1  if  musf  vanish  in 
some  domain  of  z.  Now  eilher  solution  in  Eq.  (67)  is  not  valid  near  qi(z)  0.  This  means  fhal  when  qi(z) 
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is  >  0  and  changes  to  <  0  in  Eq.  (65)  it  must  vanish,  and  by  the  preceding  (qualitative)  discussion  regarding 
the  solutions  in  Eq.  (67),  the  corresponding  qualitatively  similar  solutions  in  Eq.  (65)  shall  also  be  invalid 
near  the  zeros  defined  by  qj  (z)  =  0. 

Examining  the  solution  in  Eq.  (67)  one  notes  that  for  a  uniform  solution  valid  near  the  zeros  defined 
by  qi(z)  =  0  and  also  when  qi(z)  =  ±1,  the  functional  behavior  of  ^(z)  defined  by  Eq.  (66)  needs  fo 
be  expressible  in  terms  of  non- elementary  functions  which  would  reduce  to  the  corresponding  elementary 
functions  in  the  appropriate  regions  of  validity.  Eor  the  simplest  differential  equation  defined  in  Eq.  (66)  a 
generalized  form  would  be 

— 5 — koZ'^  =  0.  (68) 

dz 


The  preceding  is  Airy’s  differential  equation  and  can  be  solved  by  using  Eaplace  transforms  [42,  pp.  50-53]. 
The  contour  integral  representation  of  the  solution  quoted  from  Ref.  [42,  p.  52,  Eq.  2.5.9]  reads 


Bi{s)  =  j[e  ^^'^/^Ai[se  —  e^^^/^Ai{se^'^'^/^)].  (69) 


The  contour  C  in  Eq.  (69)  begins  at  o°e  and  ends  at  complex  x  plane.  The  solution  to 

Eq.  (68)  is  then  expressible  in  terms  of  Airy  functions  that  read 


fA  (z)  =  AoA/(z)  +  BoR/(z) . 


(70) 


The  various  formulas  for  numerical  calculation  of  Airy  functions  can  be  found  in  Ref.  [91].  In  Eq.  (70) 
the  constants  Ao,Bo  have  absorbed  the  wavenumber  Eq.  The  dominant  (asymptotic)  behavior  of  the  Airy 
functions  [17,  p.  373;  91]  are  as  below 


(71) 

The  preceding  asymptotic  expansions  can  be  obtained  [42,  chap.  7]  starting  with  the  contour  integral  repre¬ 
sentation  in  Eq.  (69).  Observe  that  in  light  of  Eq.  (71)  the  result  in  Eq.  (70)  can  exhibit  both  exponential  and 
oscillatory  behavior  and  hence  qualitatively  the  two  different  solutions  given  in  Eq.  (67)  appear  as  valid  in 
their  respective  regions.  However,  the  general  solution  in  Eq.  (70)  is  valid  near  the  turning  points  while  the 
various  approximations  in  Eq.  (71)  are  not.  This  suggests  that  the  one-dimensional  turning  point  problems 
for  a  continuously  stratified  media  will  contain  non-elementary  functions  as  a  part  of  its  solution.  These 


A/(z) 

Bi{z) 

A/(z) 

Bi{z) 


exp(^z3/2) 

2V^zi/4 

exp(^z^/2) 


,  when  z  ^  -|-oo, 
,  when  z  ^  -|-oo, 


1 


\/7t(-z)l/4 


sin[|(-z)^/^-hy],  when 


1 


z  — >  — 


-^=p^cos(-{-z)=/^+^],when 


Z  ^  - 


44 


Chatterjee  and  Walker 


non-elementary  functions  are  valid  across  turning  points  and  rays  in  specific  directions  that  pass  through 
these  turning  points  in  a  complex  plane  are  known  as  “Stokes  lines”.  Across  the  Stokes  lines  the  WKB 
solution  (which  contains  elementary  functions  only,  as  in  Eq.  (65))  changes  discontinuously  and  the  Phase- 
Integral  solution,  which  apparently  is  a  generalization  of  the  conventional  WKB  solution,  remains  valid. 
These  ideas  are  further  pursued  in  Refs.  [44,  88-90,  98]  and  there  exists  a  whole  new  area  of  research  in 
determination  of  numerically  accurate  solutions  valid  at  and  near  Stokes  lines.  The  results  of  this  chapter 
are  comprehensively  discussed  next. 

3.4  Summary  and  Discussion 

The  one-dimensional  Green’s  function,  as  a  solution  to  Eq.  (59),  can  be  formally  represented  as 


G^(z,z';p)  = 


Jfi(z<)^(z>) 

p{z')W{z') 


w{z!)  =  .jei{t) 


dJk2{t) 

dt 


■Mt) 


d.yk\{t) 

dt 


t=z' 


(72) 


following  the  general  form  given  in  Ref.  [38,  p.  65,  Eq.  27].  The  eigenfunctions  ^=1.2(2)  and  the  Wron- 
skian  W (z')  are  assumed  to  be  uniquely  determined.  Then,  together  with  Eqs.  (65),  (59),  the  electric  or 
magnetic  potential  Green’s  functions  in  Eq.  (57)  can  be  determined.  Note  that  the  feature  of  the  overall 
solution  given  by  Eq.  (57)  is  that  z-stratification  is  separated  from  the  large  lateral  separations  that  are  con¬ 
tained  as  arguments  of  the  Hankel  function  in  Eq.  (57).  This  feature  allows  examining  independent  spatial 
approximations  to  the  potential  Green’s  function. 

Normally,  a  Phase-Integral  solution  such  as  in  Eq.  (57),  would  need  to  be  treated  by  asymptotic  methods 
and  this  entails  evaluation  of  integrals  of  the  form  [42,  p.  252,  Eq.  7.1.1] 


I(k) 


(73) 


The  recent  trends  in  research  in  uniform  asymptotic  analysis,  reviewed  comprehensively  in  Ref.  [86],  shows 
that  resurgent  analysis  [90]  and  hyperasymptotic  techniques  [94,  95;  99,  chap.  6]  are  the  current  state-of-art 
techniques  that  need  to  be  applied  to  such  electromagnetic  problems.  The  earlier  investigations  available  in 
Refs.  [35,  chap.  4;  41-44],  which  are  currently  being  used  do  not  have  such  sophistications.  The  main  reason 
is  that  across  Stokes  lines  that  emanate  through  turning  points,  the  “compound”  asymptotic  expansion  [88] 
can  be  written  generically  as 


Compound  Asymptotic  Expansion  Dominant  Asymptotic  Expansion 

-|-  Stokes  Multiplier  x  Subdominant  Asymptotic  Expansion  (74) 


The  aspect,  useful  for  numerical  calculations,  is  that  across  the  Stokes  lines  the  subdominant  and  dominant 
expansions  interchange  roles.  The  dominant  term  ceases  to  have  any  numerical  significance  while  the  sub¬ 
dominant  term  becomes  numerically  dominant.  The  Stokes  multiplier,  which  in  most  cases  is  the  Dawson 
integral,  attains  a  maximum  value  of  1  in  the  neighborhood  of  the  Stokes  line.  As  the  observation  point 
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moves  away,  the  subdominant  term  becomes  insignificant  again,  because  the  (Stokes  multiplier)  Dawson 
integral  becomes  numerically  small.  These  changes  to  the  numerical  value  of  the  Dawson  integral  happens 
very  rapidly  at  and  around  the  Stokes  lines. 

Recent  developments  with  error  bounds  [94]  in  steepest  descent  methods  show  that  expansions  to 
Eq.  (73)  give  errors,  (f,  that  are 


j — pv+I  Poincare  type  asymptotic  expansions  [42] ; 

S’  oc  exp(— |k|A)  for  Stokes  type  asymptotic  expansions,  and 

oc  exp(— 2.386|k|A)  from  Ref.  [94].  (75) 


The  information  in  Eqs.  (74)  and  (75)  is  extremely  critical  for  numerical  work  related  to  Green’s  functions 
in  layered  media. 


4.  QUALITATIVE  REVIEW  OF  TECHNIQUES  AND  SCOPE  FOR  FURTHER  WORK 

In  this  chapter  some  of  the  state-of-art  techniques  for  radiation  of  sources  in  a  piecewise  discrete  mul¬ 
tilayer  or  uniaxially  continuously  stratified  media  terminated  in  a  PEC -backed  groundplane,  are  discussed 
qualitatively.  The  features  of  the  Sommerfeld  and  Phase-Integral  approaches  are  also  compared  with  the 
others. 

4,1  Comparative  Review  of  Some  Methodologies 

Eor  calculating  fields  in  a  layered  media,  the  integral  equation  approach  normally  takes  preference  over 
differential  equation  methods  [33,  34]  because  of  the  stability  of  the  integral  equation  operator.  However,  it 
appears  that  an  analytic  form  of  the  Green’s  function  is  still  necessary.  The  results  in  Chapter  2,  in  particular 
Eqs.  (14-16),  indicate  that  as  the  observation  point  moves  across  the  layers  the  Green’s  function  changes 
its  mathematical  form.  As  the  numbers  would  increase,  the  nature  of  the  Green’s  function  would  get  more 
complicated.  This  observation  also  indicates  that  a  straightforward  application  of  the  Einite  Difference 
Time  Domain  (EDTD)  method  becomes  almost  impossible  even  for  modest  number  of  layers  and  at  high 
frequencies.  This  is  due  to  the  fact  that  EDTD  method  typically  uses  the  discretization  of  the  layers  into 
cubic  cell  sizes  of  X  /20  side.  Smaller  size  cells  would  give  better  accuracy.  Thus,  as  the  frequency  increases 
the  number  of  cells  can  increase  to  a  number  which  would  render  computations  with  EDTD  impossible  in 
practical  terms.  Alternate  methods  such  as  the  Discrete  Complex  Image  Method  (DCIM)  [66-68],  real-axis 
integration  schemes  [56-58,  47-51],  and  asymptotic  techniques  [31,  32,  63-65,  80-85]  form  the  major  thrust 
of  the  research  work.  Other  analytical  techniques  [69-73]  have  been  gainfully  utilized,  but  these  are  not 
very  different  from  the  ones  cited  above.  In  addition,  some  of  the  recent  work  such  as  Ref.  [73]  appear  to 
solve  the  problem  in  terms  of  incomplete  Hankel-Eipschitz  integrals  which  are  quite  involved  for  numerical 
calculations.  Eurthermore,  for  the  problem  at  hand,  see  Eig.  2,  a  reliable  solution  at  extremely  large  (lOOOX) 
and  moderate  (lOX)  distances  is  required.  It  is  not  known  whether  the  techniques  in  Refs.  [69-73]  would 
be  directly  applicable,  and  if  they  would  have  any  advantage  over  the  existing  asymptotic  methods.  The 
features  of  some  of  the  above-mentioned  methods  are  discussed  below: 
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Real-axis  Integration  [56-58,  47-51]:  This  method  has  been  applied  to  the  ealeulation  of  eontaining 
Sommerfeld  integrals  for  single  and  double  layer  PEC-baeked  substrates  in  Chapter  2  (numerieal  re¬ 
sults  for  single  layer  are  ineluded).  The  most  important  aspeet  is  to  ealeulate  the  pole  loeations  on  the 
proper  Riemann  sheet  and  evaluate  the  tail  part  of  the  speetral  integral  in  elosed-form.  The  method 
ean  be  applied,  in  prineiple,  to  multiple  layers.  The  primary  diffieulties  are  that  if  the  number  of  layers 
inerease  beyond  two,  loeation  of  the  poles  and  ealeulation  of  the  residues  beeome  extremely  eumber- 
some.  However,  the  evaluation  of  the  Sommerfeld  integral  mostly  exaet  -  with  no  approximations 
to  the  Bessel  funetion,  unlike  in  Refs.  [51,  57].  Applieation  of  the  method  in  Ref.  [56]  to  eontinu- 
ously  stratified  media,  with  the  Sommerfeld  tail  part  evaluated  in  elosed-form,  appears  feasible  only 
if  sueh  a  layered  media  is  diseretized  into  a  finite  number  of  layers  with  pieeewise  eonstant  eleetrieal 
parameters. 

DCIM  [66-68]:  The  feature  of  the  DCIM  is  that  it  utilizes  the  Sommerfeld-Weyl  identity  and  ean  express 
the  Green’s  funetion  in  terms  of  some  few  exponentials  ineluding  the  image  effeets  [66,  Eqs.  12- 
14].  However,  the  number  of  surfaee  waves  poles  inereases  as  the  substrate  thiekness  inereases  and 
these  effeets  need  to  be  ineluded  [66,  Eqs.  15-16].  The  DCIM  method’s  primary  advantage  is  that 
it  expresses  the  Green’s  funetion  as  a  sum  of  finite  number  of  terms.  However,  an  examination  of 
the  reeent  work  in  Ref.  [68]  indieates  that  improvements  in  redueing  the  fill-time  of  the  matrix  in  a 
MoM  solution  have  been  aehieved  and  henee  the  new  methodology  of  the  2D-DCIM  appears  very 
eompetitive  to  solve  multilayer  problems.  It  must  be  noted  that  to  date  there  is  no  applieation  of  the 
DCIM  to  eontinuously  stratified  media.  Indeed,  eonfinuously  slralified  media  is  diserefized  info  a  sef 
of  pieeewise  eonsfanf  layers  and  fhen  DCIM  ean  be  applied  as  shown  in  Ref.  [68].  The  extension  fo 
DCIM  fo  eonfinuously  slralified  media  appears  lo  be  a  ehallenging  lask. 

Asymptotic  Techniques  [31,  32,  63-65,  80-85]:  These  are  by  far  the  most  useful  methods  beeause  they  ean 
result  in  signifieant  reduetion  of  eomputation  times  [63].  The  diffieulty  is  the  limitation  in  applieabil- 
ity.  Eor  example,  the  results  in  Ref.  [63]  are  for  a  single  layer  substrate,  but  they  are  valid  for  souree 
and  observer  loeated  at  the  air-dieleetrie  interfaee.  The  reason  for  this  limitation  is  the  nature  of  the 
asymptotie  reduetion  employed.  In  Ref.  [63]  the  modified  Sommerfeld  integral,  via  approximation  of 
the  Hankel  funetions,  are  evaluated  for  the  speeial  ease  of  “simple  pole  elose  to  a  saddle  point”.  Now 
as  the  observation  point  moves  from  inside  the  substrate  into  the  air,  many  other  analytie  singularities 
eome  into  a  eoalesenee.  The  result  is  that  sueh  a  reduetion  of  the  Sommerfeld  integral  as  in  Ref.  [63] 
becomes  invalid  for  arbitrary  source  and  observer  locations.  Similar  remarks  apply  to  Ref.  [64].  The 
situation  is  clearly  depicted  in  Refs.  [31,  32].  Here  it  is  shown  that  due  to  existence  of  branch  points 
and  poles,  which  can  come  close  to  saddle  points  the  resulting  uniform  asymptotic  evaluation  shall 
contain  more  complicated  non-elementary  functions  such  as  Weber  parabolic  cylinder  function.  Re¬ 
gardless,  application  of  some  form  of  asymptotic  techniques  appear  fruitful  when  relative  geometrical 
locations  between  source  and  observer  become  very  large. 

Phase-Integral  Methods  [16-20]  are  a  part  of  asymptotic  techniques  that  can  be  applied  to  obtain 
the  Green’s  function  for  a  continuously  stratified  media  (this  has  been  discussed  in  Chapter  3).  The 
method,  when  compared  against  a  formulations  for  discretized  layered  media  with  piecewise  con¬ 
stant  characteristics,  yields  a  single  Green’s  function  that  is  valid  for  arbitrary  observation  and  source 
locations.  However,  while  this  technique  has  been  successfully  to  problems  of  acoustic,  optical  and 
seismic  waves  propagation  in  layered  media  [2-5,  7,  28],  its  application  to  electromagnetic  wave  prop¬ 
agation  in  layered  media  with  an  emphasis  on  antenna  problems  appears  to  be  limited  [22,  24,  25].  The 
reason  can  perhaps  be  traced  back  to  the  apparent  mathematical  difficulties  in  accurately  evaluating 
integrals  across  Stokes  lines  and  caustics.  However,  with  the  advent  of  newly  developed  asymptotic 
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methods  [86-100],  the  phase-integral  solutions  to  the  layered  media  Green’s  funetion  should  not  be 
diffieult  to  ealeulate  numerieally. 

Ray-Mode  Methods  [80-85]  are  partieularly  useful  for  antenna  array  problems  where  the  sourees 
maybe  embedded  inside  layered  media  -  as  in  this  situation.  The  method  utilizes  both  mode  and  rays 
for  a  problem-matehed  Green’s  funetion  [80].  The  use  of  infinite  or  finite  Poisson  sum  formula  for 
the  array  Green’s  funetion  is  shown  to  reduee  mode  to  rays  and  viee  versa.  Eaeh  deseription,  ray  and 
mode,  is  used  in  appropriate  regions  that  finally  yields  besf  eonvergenee.  The  mefhod  apparenfly  ean 
be  applied  fo  inhomogeneous  media  [83],  and  needs  fo  be  furfher  examined  in  eonfexf  of  fhis  problem. 

4.2  Discussion  and  Scope  for  Further  Investigations 

Some  of  the  important  future  tasks  for  this  investigation  are  identified  and  itemized  below.  These  are  by 

no  means  exhaustive  but  appear  relevant  from  the  investigations. 

Task  1:  In  this  investigation,  limited  by  resourees,  a  method  to  ealeulate  Sommerfeld  integral  by  performing 
integration  along  the  real  axis  has  been  earried  out.  However,  as  shown  in  Chapter  2,  the  Bessel 
funetion  Jo(^p)  was  not  approximated  to  ealeulate  the  tail  part.  The  Sommerfeld  identity  (derived  in 
Appendix  A)  was  used.  The  method  of  integration  is  an  improvement  over  the  teehniques  in  Refs.  [51, 
57]  but  its  full  effieaey  still  needs  to  be  realized.  Beeause  it  ean  be  applied  to  diseretized,  pieeewise 
eonstant  multilayer  media  this  work  needs  to  be  earried  out.  This  effort,  in  the  least  shall  provide  a 
referenee  solution  against  whieh  results  from  many  other  teehniques  ean  be  eompared.  Extension  to 
a  two-layer  media  appears  straightforward,  though  the  pole  loeations  ean  eause  severe  problems. 

Task  2:  The  2D-DCIM  needs  to  be  applied  to  multilayered  media  and  will  also  provide  a  useful  parallel 
referenee  solution  like  the  real-axis  integration.  There  is  no  antieipated  diffieulty  beeause  it  has  been 
applied  to  two-layer  media  [68]  (as  a  separate  referenee  solution  the  EDTD,  for  limited  eases,  needs 
to  be  utilized  to  obtain  referenee  solutions  that  ean  be  eompared  against  the  2D-DCIM  and  other 
methods). 

Task  3:  Asymptotie  teehniques,  whieh  ean  inelude  straightforward  uniform  asymptotie  expansions,  phase- 
integral  methods  and  ray-mode  teehniques  appear  to  be  the  best  approaehes  for  both  multilayered 
and  eontinuously  stratified  media.  In  applying  the  asymptotie  teehniques,  eare  must  be  exereised  to 
ineorporate  reeent  researeh  in  this  area  sueh  as  aeeurate  numerieal  evaluation  of  the  Green’s  funetion 
aeross  Stokes  lines. 

Task  4:  Beeause  asymptotie  teehniques,  DCIM,  and  real-axis  integration  apply  to  diseretized  layered  me¬ 
dia,  validation  of  the  various  teehniques  should  eonsider  the  diseretized  layered  media  as  an  important 
benehmark  problem.  Next,  the  phase-integral  method  should  be  formulated  for  eontinuously  stratified 
ease  and  eompared  against  the  diseretized  ease.  Similar  remarks  apply  for  the  ray-mode  teehnique. 
As  a  eonsequenee,  a  very  eritieal  task  is  to  aseertain  with  some  degree  of  aeeuraey  the  effeets  of 
rapidly  varying  permittivity,  e(z),  and  permeability,  |4(z),  profiles  when  the  eontinuous  variation  is 
diseretized. 

Task  5:  In  this  report  only  the  eomponent  was  studied  for  a  horizontal  eleetrie  eurrent  element.  The 
results  for  an  arbitrarily  oriented  Hertzian  eurrent  element  needs  to  be  found  to  eomplete  the  entire 
problem.  All  the  preeeding  methods  needs  to  be  applied  to  obtain  solutions  for  the  arbitrarily  oriented 
Hertzian  dipole  ease. 
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In  addition,  the  methods  should  also  be  applied  to  eonformal  arrays  embedded  inside  sueh  layered  media. 
The  survey  of  literature  [9,  eh.  4;  10]  in  this  area  indieates  that  there  exists  a  paueity  of  useful  numerieal 
modeling  teehniques  for  sueh  elass  of  array  radiators. 


5.  CONCLUSION 

In  this  report  radiation  of  sourees  inside  layered  media,  baeked  by  a  perfeet  eleetrie  eonduetor  were 
studied.  This  problem  has  praetieal  applieations  to  antennas  embedded  inside  eomposite  material  media  for 
many  applieations.  The  fields  due  to  a  horizontal  eleetrie  eurrent  element  were  studied  by  ealeulating  the 
single  Gjx  eomponent  of  the  dyadie  Green’s  funetion  due  to  a  horizontal  eleetrie  eurrent  element.  The  tra¬ 
ditional  real-axis  integration  method  was  employed  with  a  new  formulation  for  ealeulating  the  Sommerfeld 
integral  tail.  This  new  method  eontains  effects  of  the  intrinsic  wavenumber  of  the  layer  and  hence  appears 
physically  more  suitable  for  understanding  effects  of  the  various  layers  when  the  observation  point  moves  in 
a  vertical  direction  from  one  layer  to  the  other.  In  addition,  a  new  algorithm  was  developed  to  calculate  the 
proper  surface  wave  poles  for  electrically  thick  substrates.  The  discontinuities  in  the  integrand  of  the  Som¬ 
merfeld  integral  were  eliminated  upon  subtraction  of  the  residues  at  these  poles,  and  subsequent  numerical 
integration  posed  no  difficulties.  The  calculations  were  also  performed  via  the  computer  software  MATH- 
EMATICA  through  direct  numerical  integration  without  any  preceding  modifications.  The  results  from  the 
two  approaches  showed  that  the  real-axis  integration  algorithm  was  at  least  several  times  faster  compared 
to  direct  numerical  integration  via  MATHEMATICA  when  the  lateral  separation  between  the  source  and 
observer  points  increased,  with  no  loss  of  accuracy.  Additionally  the  problem  of  stratified  media  was  sfud- 
ied  via  fhe  Phase-Infegral  mefhod.  If  has  been  shown  fhaf  fhe  Phase-Infegral  mefhod,  wifh  ifs  affendanf 
asympfofic  expansions,  has  fhe  pofenfial  for  providing  much  improved  solufions  fo  fhe  Green’s  functions  for 
confinuous  media  by  including  effecls  across  fhe  Sfokes  lines.  Einally,  some  recommendafions  for  furfher 
work,  based  on  fhe  presenf  invesfigafions,  have  been  included. 
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CONVERSION  OF  TWO-DIMENSIONAL  SPECTRAL  INTEGRALS  TO 

SOMMERFELD  FORMS 

Shaun  Walker 


In  this  appendix,  the  eonversion  of  the  two-dimensional  plane-wave  speetral  (PWS)  integral  from  Ref.  [Al, 
p.  188,  Eq.  4-124]  is  eonverted  to  the  well-known  Sommerfeld  identity.  Following  the  Sommerfeld  identity, 
some  algebraie  manipulations  are  further  performed  to  obtain  the  desired  relationships  for  evaluating  the 
“tail”  part  of  the  Sommerfeld  integral  in  exaet,  closed  form.  These  details  are  included  in  this  appendix.  The 
PWS  spectral  integral  reads 
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The  following  polar  transformations  on  the  spectral  and  spatial  variables  are  utilized:  kx  =  '^  cos  a ,  ky  = 
^  sin  a  and  x  =  p  cos(|),  y  =  p  sin(|).  This  transformation  yields  dkxdky  =  ^d^da.  The  region  defined  by 
— oo  <  {kx,ky)  <  -|-oo  is  mapped  onto  the  region  0  <  a  <  271  and  0  <  ^  <  -|-oo.  Consequently  Eq.  (Al) 
reduces  to  the  form 
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In  Eq.  (A2)  the  following  relationship  for  the  generating  function  for  Bessel  function 
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is  used  to  perform  the  following  reduction 
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Since 
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j+im(a-(|))jp^  _ 


0,  when  m^O, 
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it  follows,  using  Eqs.  (A5)  and  (A4)  into  Eq.  (A2)  that 
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Next,  it  follows  readily  from  Eqs.  (A6)  and  (Al)  the  well-known  Sommerfeld  identity 
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When^  >  k,  then  — y'V^^^^-^inEq.  (A7).  In  both  Eqs.  (Al)  and  (A7),  r=y^x^  +  y2  -|-z^  = 

+  z2.  With  the  above  changes  in  Eq.  (A7)  and  alternate  form  of  the  Sommerfeld  identity  reads 
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Eq.  (A8)  is  known  as  the  classical  form.  In  Eqs.  (Al),  (A7),  and  (A8)  it  is  assumed  that  z  >  0  for  the 
integrals  to  converge.  Next  we  derive  the  identities/relationships  used  for  calculating  the  “tail”  part  of  the 
Sommerfeld  integral  along  real  axis. 

Differentiating  both  sides  of  Eq.  (A7)  w.r.t  p ,  and  using  the  relationships 
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In  Eq.  (A9),  it  maybe  recalled  that  r  =  y^p^  +  z^.  Next,  we  differentiate  Eq.  (A9)  w.r.t  z,  and  substitute  the 
following  relationship 
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to  obtain  the  following 


-ykr 


-[3  +  37kr-(kr)2]  =  <^ 


/  k  when^  <  k, 

b 

■h®®  / - o 

/  when^>k 


(AlO) 


A  slightly  different  variant  of  Eq.  (AlO)  was  used  in  Ref.  [A2,  Eq.  A4].  This  form  can  also  be  obtained  from 
Ref.  [A3,  p.  694,  Eq.  6.623.2]  and  reads 


+  00 

y"^2ji(^p)exp(-^z)d^=^.  (All) 

0 

This  form  can  be  used  to  calculate  a  different,  final  asymptotic  form  for  the  slab  function  )  given 

explicitly  in  Eq.  (13).  The  form  in  Eq.  (All)  can  be  obtained  as  a  special  case  from  Eq.  (AlO)  by  setting 
k  =  0  therein  for  ^  >  k  case.  The  accuracy  of  using  both  forms  in  numerical  calculation  of  Sommerfeld 
integral  tails,  however,  remains  an  open  question. 
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Appendix  B 


ADDITIONAL  RESULTS  FOR  THE  INTEGRAND  IN  THE  Gzx  GREEN’S  FUNCTION 
FOR  SINGLE  LAYER  PEC-BACKED  SUBSTRATE 


Fig.  B1  —  Behavior  of  integrand  in  Eq.  (3)  for  observation  point  in  air  (region  #  0)  in  Fig.  3  for 
p  =  1000?i,z  =  +|,tan5  =  0.0001,  d  =  A,  and  |eH  |  =  9.8 
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Fig.  B2  —  Behavior  of  integrand  in  Eq.  (3)  for  observation  point  in  air  (region  #  0)  in  Fig.  3  for 
p  =  lOOOX,  z  = +^,  tan5  =  0.0001,  d  =  A,  andle^ll  =9.8 


Fig.  B3  —  Behavior  of  integrand  in  Eq.  (3)  for  observation  point  in  air  (region  #  0)  in  Fig.  3.  All  data 
same  as  in  Fig.  B2,  except  for  the  lossless  case  tan  8  =  0. 


Fig.  B4  —  Behavior  of  integrand  in  Eq.  (4)  for  observation  point  in  substrate  (region  #  1)  in  Fig.  3 
for  p  =  lOOOA,,  z  =  —  tanS  =  0.0001,  d  =  A,  and  |e^i|  =  9.8 


Fig.  B5  —  Behavior  of  integrand  in  Eq.  (4)  for  observation  point  in  substrate  (region  #  1)  in  Fig.  3 
for  p  =  lOOOA,,  z  =  —  tan5  =  0.0001,  d  =  >,  and  |eri|  =  9.8. 
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R'tyi(f)] 


Fig.  B6  —  Behavior  of  the  integrand  in  the  Sommerfeld  integral  for  the  Gj;t(p ,  z,  (|) )  Green’s  function 
in  Eqs.  (32)  and  (33).  Here  d  =  ^,  p  =  lOOOA,,  z  =  —j,  tan6  =  0,  and  le^i  |  =  9.8.  The  observation 
point  is  inside  the  substrate.  The  imaginary  part  of  the  integrand  identically  vanishes. 


Fig.  B7  —  Behavior  of  the  integrand  in  the  Sommerfeld  integral  for  the  G2;c(p ,  z,  () )  Green’s  function 
in  Eqs.  (32)  and  (33).  Here  d  =  A,,  p  =  lOOOA,  z  =  +^,  tanS  =  0,  and,|eri|  =  9.8.  The  observation 
point  is  outside  the  substrate.  The  real  part  of  the  integrand  identically  vanishes. 


Appendix  C 


MATHEMATICA  SCRIPTFILE  (CODE  LISTING)  FOR  PROPER  SURFACE  TM  WAVE 

POLE  LOCATIONS 

Shaun  Walker 


(* 

List  of  physcal  constants 


cO  -  Speed  of  light  in  vacuum 

f  -  operating  frequency 

lam  -  wavelength 

kO  -  free-space  wavenumber 

tand  -  dielectric  loss 

er  -  relative  permittivity  of  slab 

d  -  thickness  of  slab 

kl  -  wavenumber  in  slab 


OtherVariables 


numZEROS  -  number  of  zero  locations  of  the  denominator 

in  the  integrand  (i.e.  number  of  zeros  of  DTM) 
orderOfTS  -  order  of  the  Taylor  series  used 
correctZEROS  -  this  is  used  inside  of  a  loop  to  count  the 
number  of  correct  zeros  found  by  the  algorithm,  when  variable 
is  equal  to  numZEROS  then  this  means  that  all  roots  where  found 
zeroLIST  -  list  of  zeros  that  contains  the  zeros  of  the 
Taylor  series  that  correspond  to  the  proper  surface  wave  poles 
testROOTS  -  a  list  of  length  orderOfTS,  it  contains  roots 
taken  from  the  zeros  of  the  truncated  Taylor  polynomial  which 
include  the  effects  of  all  modes,  these  roots  will  be  tested  to  see 
which  ones  correspond  to  the  surface  wave  poles 


Built-in  Mathematica  subroutines 


Sqrt [ ] 
Floor [ ] 
Re[] 
lm[] 

Abs  [  ] 
NRoot [ ] 


-  takes  the  square  root  of  its  argument 

-  truncates  the  decimal  portion  of  a  real  number 

-  takes  the  real  part  of  a  complex  number 

-  takes  the  imaginary  part  of  a  complex  number 

-  takes  the  absolute  values  of  a  complex  number 

-  numerically  finds  roots  of  a  polynomial 
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FindRoot [ ]  -  numerically  finds  root  of  a  user  defined  function 


Mathematica  Constants 


Pi  -  pi=3 . 1415459045 
I  -  square  root  of  -1 
* ) 

{*  Subroutine  used  to  evaluate  the  mode  function  *) 

DTM[d_,  er_,  k0_,  v_]  :=  Module [{D}, 

D  =  er*v*Cos  [ Sqrt  [  (d*k0 ) ''2*  (er-1 ) -v''2 ]  ] 

-  Sqrt  [  (d*k0)  "2*  (er-l)-v''2]*Sin[  Sqrt  [  {d*k0 )  "  2  *  (er-1 )  -  v"  2  ]  ]  ; 

D] 

{*  define  constants  *) 

cO  =  2 . 997924*10''8; 

f  =  5*10''9; 

lam  =  cO/f; 

kO  =2  Pi/lam; 

tand  =  0.0001; 

er  =9.8  { 1  -  I  tand) ; 

d  =  lam; 

kl  =  Sqrt [er]  kO; 

numZEROS  =  Floor [d  kO  Sqrt [Re [er] -1] /Pi] +1; 

orderOfTS  =  0; 
correctZEROS  =  0; 

{* 

while  the  number  of  correct  zeros  found  by  the  algorithm  is  less 
than  the  number  of  zeros  of  DTM  due  to  proper  surface  wave  poles, 
then  keep  looking 
* ) 

While [ correctZEROS  <  numZEROS, 

{*  reset  the  number  of  correct  zeros  found  &  list  of  zeros  *) 
correctZEROS  =  0; 

zeroLIST  =  Table [0,  [i,  numZEROS}]; 

{* 

increase  the  order  of  the  Taylor  series,  and  insure  that  the  order 
of  the  Taylor  series  taken  will  be  greater  than  the  number  of 
surface  wave  poles  that  need  to  be  found 
* ) 

If [ (orderOfTS++)  <  numZEROS,  orderOfTS  =  numZEROS]; 

{*  get  Taylor  series  of  DTM  and  find  all  zeros  of  the  series  *) 
testROOTS  =  NRoots [Normal [Series [DTM[d,er,k0,v] , {v, 0, orderOfTS } ] ] 

==0, v] ; 

{*  for  each  of  the  test  roots  found  from  the  truncated  Taylor  series  *) 
For[q  =1,  q  <=  Length [testROOTS ] ,  q++, 

{*  gets  the  next  test  root  *) 

If [orderOfTS  >1,  vr  =  testROOTS [ [q] ] [ [2] ] ,  vr  =  testROOTS [ [1] ] ] ; 
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{*  evaluate  DTM  for  test  root*) 
test  =  DTM[d,  er,  kO,  vr] ; 

{*  is  the  test  close  enough  to  zero  and  is  this  a  proper  SWP  ?  *) 

If [Abs [test ] <=1  &&  Re[vr]>=0  &&  Re [vr] <=d  kO  Sqrt [Re [er ] -1 ]  && 

Im [ vr ] <=0 , 

{*  then  increment  the  number  of  correct  zeros  that  have  been  found  *) 
correct ZEROS ++; 

{*  add  the  test  root  to  the  list  of  correct  zeros  found  *) 
zeroLIST [ [correctZEROS] ]  =  vr 

] ; 


{*  use  the  inital  guess  for  each  root  to  find  the  exact  root  *) 
exact ZEROS  =  FindRoot [DTM [d,er,kO,v] , {v, zeroLIST}, AccuracyGoal->IO ] ; 
exactZEROS  =  Sqrt [ (exactZEROS [ [I] ] [ [2] ] / (d  kO) ) "2  +  I] ; 

Print [ "Number  of  zeros  is  ",  numZEROS, 

Print ["Order  of  the  Taylor  series  taken  is  ",  orderOfTS,  " . " ] ; 

Print ["Zero  locations  are:  ",  exactZEROS]; 
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Appendix  D 


MATHEMATICA  SCRIPTFILE  (CODE  LISTING)  FOR  CALCULATING  THE  Gzx 
GREEN’S  FUNCTION  FOR  PEC-BACKED  SINGLE  LAYER  MEDIA 

Shaun  Walker 


{* 

There  are  6  subroutines,  two  of  these  functions  are  for  the 
evaluation  of  the  residues,  and  the  other  4  equations  are 
for  the  calculation  of  Gzx. 

RpO  -  Calculates  the  pth  residue  in  region  0 
Rpl  -  Calculates  the  pth  residue  in  region  1 
zLessThanZero  -  Computation  of  Gzx  for  z<0,  in  Region  0 
zEqualZeroRl  -  Computation  of  Gzx  for  z=0  using  region  1  expression 

zEqualZeroRO  -  Computation  of  Gzx  for  z=0  using  region  0  expression 

zGreaterThanZero  -  Computation  of  Gzx  for  z>0,  in  Region  1 

* ) 

RpO [Rho_, z_, d_, er_, k0_, Xi_]  :=  Module [ {R} ,  kl=Sqrt [er ] *k0; 

Num  =  Xi'Z  *  BesselJ [ 1 , Rho*Xi ]  *  Sqrt [kl ' 2-Xi "2 ]  * 

Sin  [d*  ( Sqrt  [kl''2-Xi]''2])]*  E"  (-z*Sqrt  [Xi''2-k0''2]  )  ; 

Dp  =  -I*Xi*Cos [d*Sqrt [kl "E-Xi ' 2 ] ] * (d+er/Sqrt [Xi " 2-kO "2 ] ) 

-  I*Xi*Sin  [d*Sqrt  [kl''2-Xi''2]  ]  *  (l+d*er*Sqrt  [Xi''2-k0''2]  ) 

/ Sqrt [kl " 2-Xi " 2 ] ; 

R=Num/Dp ; 

R] 

Rpl [Rho_, z_, d_, er_, k0_, Xi_]  :=  Module [{R},  kl=Sqrt [er ] *k0; 

Num  =  Xi  "  2  *BesselJ  [  1 ,  Rho*Xi  ]  *  (-1  *Sqrt  [Xi ''2-kO  "  2  ]  ) 

*  Cos [ (d-Abs [ z ] ) *Sqrt [kl " 2-Xi ' 2 ] ] ; 

Dp  =  -I*Xi*Cos  [d*Sqrt  [kl ''2-Xi''2  ]  ]  *  (d+er / Sqrt  [Xi  " 2-kO  " 2  ]  ) 

-  I*Xi*Sin [d*Sqrt [kl " 2-Xi "2 ] ] /Sqrt [kl " 2-Xi "2 ]  * 

( l+d*er*Sqrt  [Xi''2-k0''2]  )  ; 

R=Num/Dp ; 

R] 

zLessThanZero [Rho_, z_, d_, erl_, k0_, XiA_]  ;=  Module [{ Integral } , 
kl  =  Sqrt [erl ] *k0; 
r  =  Sqrt  [Rho''2  +  z''2]  ; 

lOtoKO  =  NIntegrate [Xi ' 2 *BesselJ [ 1 , Rho*Xi ] *Sqrt [kO " 2-Xi ' 2 ] 

*  Cos  [  (d-Abs  [  z ]  )  *Sqrt  [kl ''2-Xi ' 2  ]  ]  /  {erl*Sqrt  [kO''2-Xi''2] 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  +  I*Sqrt  [kl''2-Xi''2]  *Sin  [d*Sqrt  [kl''2-Xi''2]  ]  ) 

-  Exp  [  -  I*Abs  [  z  ]  *Sqrt  [kl''2-Xi''2]  ]  /  (erl  +  l)  )  , 

{Xi, 0, kO } , 
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AccuracyGoal->6, 

MaxRecursion->1000, 

Method-> { GlobalAdaptive, 

MaxErrorIncreases ->200000 }  ]  ; 

IKOtoKl  =  NIntegrate[Xi''2*BesselJ[l,  Rho*Xi]  *  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [  (d-Abs  [  z ]  )  *Sqrt  [kl ''2-Xi ' 2  ]  ]  /  (erl*  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  +  I*Sqrt  [kl''2-Xi"2]  *Sin  [d*Sqrt  [kl''2-Xi''2]  )  ] 

-  Exp  [  -  I*Abs  [  z  ]  *Sqrt  [kl''2-Xi''2]  ]  /  (erl  +  1)  ) 

-  Sum (  Rpl [Rho, z, d, erl, kO, Xip [ [p] ] ]  /  (Xi-Xip[[p]]  ,  {p=l,P}), 

{Xi,  kO, Re [kl]  } , 

AccuracyGoal->6, 

MaxRecursion->1000, 

Method-> { GlobalAdaptive, 

MaxErrorIncreases ->200000 } ] 

+  Sum[  Rpl [Rho, z, d, erl, kO, Xip [ [p] ] ] 

*  Log [ (kl-Xip [ [p] ] ) / (Xip [ [p] ] -kO) ] -I*Pi  ,  {p=l,P}]; 

IKltoXiA  =  NIntegrate [Xi " 2*BesselJ [ 1 , Rho*Xi ]  *  (-1 *Sqrt [Xi " 2-kO "2 ] ) 

*  Cos [ (d-Abs [z] ) * (-I*Sqrt [Xi " 2-kl " 2 ] ) ] 

/  (  erl*  {-I*Sqrt  [Xi''2-k0''2]  )  *  Cos  [d*  {-I*Sqrt  [Xi''2-kl''2]  )  ] 

+  I  *  (-1  *Sqrt  [Xi  "  2-kl  "2  ]  )  *Sin  [d*  (-I*Sqrt  [Xi  ■'2-kl "  2  ]  )  ]  ) 

-  Exp  [ -Abs  [ z ]  *Sqrt  [Xi''2-kl''2]  ]  /  (erl  +  1)  , 

{Xi,Re[kl] ,XiA}, 

AccuracyGoal->6, 

MaxRecursion->1000, 

Method-> { GlobalAdaptive, MaxErrorIncreases->2 00 0  00  }  ] ; 

Itail  =  l/(erl+l)  *  Rho*z  * (3+I*3*kl*r- (kl*r) "2) 

*  Exp [-I*kl*r] /r'S; 

Integral  =  lOtoKO  +  IKOtoKl  +  IKltoXiA  +  Itail; 

Integral ] 

zEqualZeroRl [Rho_, z_, d_, erl_, k0_, XiA_]  :=  Module [{ Integral } , 
kl=Sqrt [erl] kO; 

lOtoKO  =  NIntegrate  [  (Xi''2  *  BesselJ  [  1 ,  Rho*Xi  ]  *  Sqrt  [k0''2-Xi''2] 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  /  ( erl  *  Sqrt  [  kO  "  2-Xi  '  2  ] 

*Cos [d*Sqrt [kl " 2-Xi "2 ] ] 

+  I *Sqrt  [kl''2-Xi''2]  *Sin  [d*Sqrt  [kl''2-Xi''2]  ]  , 

{Xi, 0, kO } , 

AccuracyGoal->6, 

MaxRecursion->1000, 

Method-> { GlobalAdaptive, 

MaxErrorIncreases ->200000  }  ] ; 

IKOtoKl  =  NIntegrate  [Xi''2  *  BesselJ  [  I ,  Rho*Xi  ]  *  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  /  (erl*  (-1  *Sqrt  [Xi " 2-kO ' 2  ]  ) 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  +I*Sqrt  [kl ''2-Xi''2  ]  *Sin  [d*  Sqrt  [kl " 2-Xi "2  ]  )  ] 

-  Sum[  Rpl [Rho, z, d, erl, kO, Xip [ [p] ] ]  /  (Xi-Xip [ [p] ] ) ,  {p=l,P}], 

{Xi, kO, Re  [kl]  } , 

AccuracyGoal->6, 
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MaxRecursion->1000, 

Method-> { GlobalAdaptive, MaxErrorIncreases->200000 } ] 

+  Sum [  Rpl [Rho, z , d, erl , kO , Xip [ [p] ] ] *Log [ (kl-Xip[ [p] ] ) / (Xip[ [p] ]-k0) ] 
-  I*Pi,  {p=l,P}  ] ; 

IKltoXiA  =  NIntegrate  [  (Xi ' 2  *  BesselJ  [1,  Rho*Xi]  *  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*  {-I*Sqrt  [Xi''2-kl''2]  )  ]  /  (erl*  (-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*  {-I*Sqrt  [Xi''2-kl''2]  )  ]+I*  (-1  *Sqrt  [Xi''2-kl''2]  ) 

*  Sin[d*{-I  Sqrt  [Xi  ■'2-kl  "2  ]  )  ]  )  , 

{Xi,Re[kl] ,XiA}, 

AccuracyCoal->6, 

MaxRecursion->1000, 

Method-> { ClobalAdaptive, 

MaxErrorIncreases ->200000  }  ] ; 

Itail  =  -l/{erl  +  l)  *  (XiA''2*BesselJ  [2,  XiA*Rho]  )  /  Rho; 

Integral  =  lOtoKO  +  IKOtoKl  +  IKltoXiA  +  Itail; 

Integral ] 

zEqualZeroRO [Rho_, z_, d_, erl_, k0_, XiA_]  :=  Module [{ Integral } , 
kl  =  Sqrt [erl ] kO ; 

lOtoKO  =  NIntegrate [  (k0*Sech  [u] ) "2  *  Bessel J [ 1 , Rho* (k0*Sech [u] ) ] 

*  Sqrt  [kl  "B- (kO  Sech[u])''2]  *  Sin  [d*Sqrt  [kl''2- (k0*Sech  [u]  )  "2]  ] 

*  (-k0*Sech [u] *Tanh [u] )  /  {  erl*Sqrt [kO " 2- (k0*Sech [u] ) ' 2 ] 

*  Cos  [d*Sqrt  [kl "  2- {k0*Sech  [u]  )  "2  ]  ]  +  I*Sqrt  [kl ''2- (kO  Sech[u])''2] 

*  Sin [d*Sqrt [kl " 2- {k0*Sech [u] ) "2 ] ]  ), 

[u, (kO  Sech[0]),(k0  Sech[k0])}, 

AccuracyCoal->6, 

MaxRecursion->1000, 

Method-> { ClobalAdaptive, 

MaxErrorIncreases ->200000 }  ]  ; 

IKOtoKl  =  NIntegrate  [  (Xi  "  2  *  BesselJ  [  1 ,  Rho*Xi  ]  *  Sqrt  [kl''2-Xi"2] 

*  Sin  [d*Sqrt  [kl''2-Xi''2]  ]  /  (erl*  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*Sqrt  [kl''2-Xi''2]  ]  +  I*Sqrt  [kl''2-Xi''2]  *Sin  [d*Sqrt  [kl''2-Xi''2]  ]  ) 

-  Sum[  RpO [Rho, z , d, erl , kO , Xip [ [p] ] ] / (Xi-Xip [ [p] ] ) , {p=l , P } ] , 

[Xi, kO, Re  [kl]  } , 

AccuracyCoal->6, 

MaxRecursion->1000, 

Method-> { ClobalAdaptive, MaxErrorIncreases->2 00 000 } ] 

+  Sum [RpO [Rho, z,d,erl,k0,Xip[ [p] ] ] *Log [  (kl-Xip  [  [p] ] ) /  (Xip  [  [p] ] -kO ) ] 
-I*Pi,  {p=l,P}  ]; 

IKltoXiA  =  NIntegrate[(Xi''2*BesselJ[l,  Rho*Xi ]  *  ( -I* Sqrt  [Xi''2-kl''2]  ) 

*  Sin  [d*  (-I*Sqrt  [Xi''2-kl''2]  )  ]  /  (erl  ( -I* Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*  (-I*Sqrt  [Xi''2-kl''2]  )  ]+I*  (-1  *  Sqrt  [Xi''2-kl''2]  ) 

*  Sin[d*{-I  Sqrt  [Xi  ■'2-kl  "2  ]  )  ]  )  , 

{Xi,Re[kl] ,XiA}, 

AccuracyCoal->6, 

MaxRecursion->1000, 

Method-> { ClobalAdaptive, MaxErrorIncreases->2 00 0  00  }  ] ; 
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Itail  =  -  (-1/ (erl  +  1 )  )  *  (XiA''2*BesselJ  [2,  XiA*Rho]  )  /  Rho; 

Integral  =  lOtoKO  +  IKOtoKl  +  IKltoXiA  +  Itail; 

Integral ] 

zGreaterThanZero [Rho_, z_, d_, erl_, k0_, XiA_]  : =  Module [ { Integral } , 
kl  =  Sqrt [erl ] kO ; 
r  =  Sqrt  [Rho''2  +  z''2]  ; 

1 0toK0=N Integrate  [Xi''2*BesselJ[I,Rho*Xi]  *Exp  [-I*z *Sqrt  [kO''2-Xi''2]  ] 

*  ( Sqrt  [kI''2-Xi''2]  *Sin  [d*Sqrt  [kI''2-Xi''2]  ]  /  (erl  *  Sqrt  [kO''2-Xi''2] 

*  Cos  [d*Sqrt  [kI''2-Xi''2]  ]  +  I*Sqrt  [kI''2-Xi''2]  *Sin  [d*Sqrt  [kI''2-Xi''2]  ]  ) 

+  I/(erI+I)),{Xi,0,k0} , AccuracyGoal->6, 

MaxRecursion->I000, 

Method-> { GlobalAdaptive, 

MaxErrorIncreases ->200000 }  ]  ; 

IKOtoKI  =  NIntegrate  [Xi''2  *BesselJ  [  I ,  Rho*Xi  ]  *Exp  [ -z*Sqrt  [Xi " 2-kO ''2  ]  ] 

*  (  Sqrt  [kl''2-Xi''2]  *Sin  [d*Sqrt  [kI''2-Xi''2]  ]  /  (erl*  {-I*Sqrt  [Xi''2-k0''2]  ) 

*  Cos  [d*Sqrt  [kI''2-Xi''2]  ]  +  I*Sqrt  [kI''2-Xi"2]  *Sin  [d*Sqrt  [kI''2-Xi"2]  ]  ) 

+  I/{erI+I)  )-  Sum [RpO [Rho, z , d, erl , kO , Xip [ [p] ] ] / (Xi-Xip [ [p] ] ) , {p=I , P } ] , 
{Xi, kO, Re  [kl]  } , 

AccuracyGoal->6, 

MaxRecursion->I000, 

Method-> { GlobalAdaptive, MaxErrorIncreases->200000 } ] 

+  Sum [RpO [Rho, z,d,erl,k0,Xip[ [p] ] ] * (Log [ (kl-Xip [ [p] ] ) / 

(Xip [ [p] ] -kO) ] -I*Pi) , {p=l, P } ] ; 

IKltoXiA  =  NIntegrate[Xi''2*BesselJ[l,Rho*Xi]*Exp  [-z*Sqrt  [Xi''2-k0''2]  ) 

*  ( -I* Sqrt  [Xi''2-kl"2]  *Sin[d*  ( -I* Sqrt  [Xi''2-kl''2]  )  ] 

/  (  erl  *  (-1  *  Sqrt  [Xi"2-k0''2]  )  *Cos  [d*  ( -I*  Sqrt  [Xi''2-kl''2]  )  ] 

+  I*  {-I*Sqrt  [Xi''2-kl''2]  )  *Sin[d*  (-1  *  Sqrt  [Xi''2-kl''2]  )  ]  ) 

+  1/ (erl+l) ) , 

{Xi,Re[kl] ,XiA}, 

AccuracyGoal->6, 

MaxRecursion->1000, 

Method-> { GlobalAdaptive, 

MaxErrorIncreases ->2000 00  } ]  ; 

Itail  =  -(1/ (erl+1)) *Rho*z* (3+I*3*kO*r- (k0*r) "2) *Exp [-I*k0*r] /r'S; 
Integral  =  lOtoKO  +  IKOtoKl  +  IKltoXiA  +  Itail; 

Integral ] 

Mu0  =  1.256637*10''-6; 

Epsilon0  =  8.854187*10''-12; 
c0=2 . 997924*10''8; 

f=5*10''9;  (*  operating  frequency  *) 

Lambda=c0/f;  (*  wavelength  *) 

kO=2*Pi/Lambda;  (*  free-space  wavenumber  *) 

tanDelta=0 . 0001;  (*  dielectric  loss  of  slab  *) 

erl=9 . 8* {l-I*tanDelta) ;  (*  relative  perm,  in  region  1  *) 

kl=k0*Sqrt [erl ] ;  {*  wavenumber  in  region  1  *) 

d=Lambda;  {*  slab  thicknes  *) 
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XiA=Re[kl];  (*  breakpoint  *) 

Rho  =  100*Lambda; 

Phi  =  0; 

(*  zeros  for  tand=0 . 0001,  erl=9.8,  and  d=l*Lambda  *) 

Xip  =  k0*{l. 560465260348978-1*0. 00028915689324715506, 
2.1955007373508435-1*0.00021901803692256456, 

2 . 6035344541354384-1*0 . 00018677914895957531, 

2 . 8733530563113567-1*0 . 00016999491295840665, 

3.040363182302955-  1*0.0001610010843909009, 
3.1206058917706754-1*0.00015700389753638654 
}; 

{*number  of  SWP  is  equal  to  the  length  of  the  list  of 
zeros  of  Xip  *) 

P  =  Length [Xip] ; 

{*list  of  z  values  used  to  compute  Gzx  *) 

zLIST={-l,-.95,-0.9,-.85,-.8,-.75,-.7,-.65,-.6,-.55,-.5,-.45,-.4, 
-.35, -.3, -.25, -.2, -.15,-.!, -.08, -.06, -.04, -.02, -.01,  0,. 01,. 02,. 04, 
.06, .08,.!, .15, .2, .25, .3, .4, .5, .6, .7, .8, .9, 1,1. 25, 1.5, 1.75, 2, 2. 5, 
3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10}  *  Lambda; 

{*  number  of  times  that  Gzx  needs  to  be  computed  *) 
zLENGTH=Length [zLIST] ; 

{*  initialize  Gzx  to  all  zeros  *) 

GzxA=Table [0, {i, zLENGTH+1} ] ; 

{*  for  each  values  of  z  *) 

For [m=l,m<=zLENGTH,m++,  Print ["Loop  ",m, "  of  ",zLENGTH]; 
z  =  zLIST[[m]];  (*  get  next  value  of  z  *) 

If  [z<0, 

GzxA[[m]]  =  (-I*erl)  /  (2*Pi*k0''2)  *  Cos[Phi] 

*  zLessThanZero [Rho, z, d, erl, kO, XiA] ; 

]  ; 

If [z==0, 

GzxA[[m]]  =  (-1  erl)  /  (2*Pi*kO''2)  *  Cos[Phi] 

*  zEqualZeroRl [Rho, z,d,erl,k0,XiA]  ; 

GzxA[[m+l]]  =  1/ (2*Pi*kO "E )  *  Cos [Phi] 

*  zEqualZeroRO [Rho, z,d,erl,k0,XiA] ; 

] ; 

If  [z>0, 

GzxA[[m+l]]  =  1/ (2*Pi*kO "E )  *  Cos [Phi] 

*  zGreaterThanZero [Rho, z, d, erl, kO, XiA] ; 

] 


]  / /Timing 


